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ANALYSIS  OF  THE  FINITE  ELEMENT 
RADIATION  MODEL  CODE  FOR  ANTENNA 

MODELING  * 


Marcos  C.  Medina,  Constantine  A.  Balanis,  Jian  Peng  and  Panayiotis  A.  Tirkas 

Department  of  Electrical  Engineering 
Telecommunications  Research  Center 
Arizona  State  University 
Tempe,  Arizona  85287-7206 


Abstract 

The  Finite  Element  Radiation  Model  (FERM) 
code  is  applied  to  examine  radiation  patterns 
and  input  impedance  of  various  antenna  con¬ 
figurations  on  complex  helicopters.  This  pa¬ 
per  addresses  the  strengths  and  weaknesses  of 
the  FERM  code  when  appUed  in  antenna  anal¬ 
ysis  and  is  compared  with  other  national  codes, 
such  as  the  Numerical  Electromagnetics  Code 
(NEC).  In  addition,  guidehnes  are  outlined  for 
various  applications  of  the  code. 

I.  INTRODUCTION 

The  Finite  Element  Radiation  Model  (FERM) 
code  is  based  on  the  Method  of  Moments 
and  uses  triangular  surfaces  patches  to  model 
the  geometry.  The  triangular  patches  allow 
for  enormous  flexibility  in  modehng  complex, 
curved  three-dimensional  geometries,  such  as 
hehcopters.  The  basis  functions  used  by  the 
code  are  in  accordance  with  the  Rao  et  al.  [1] 
basis  functions.  The  FERM  code  was  devel¬ 
oped  in  1987  at  the  MIT  Lincoln  Laboratory  [2], 
and  is  a  general  purpose  code  with  apphcations 
in  radar  cross  section  prediction  and  antenna 
modeling.  A  similar  code  is  the  Electromagnetic 
Surface  Patch  (ESP)  code  [3],  which  is  based  on 

‘This  work  was  supported  by  the  Advanced  Heli¬ 
copter  Electromagnetics  (AHE)  Program  and  NASA  un- 
der  Grant  NAG-1~1082. 


quadrilateral  patches.  Unlike  the  FERM  and 
ESP  codes,  the  NEC  is  based  on  a  wire  grid 
model  of  the  geometry.  The  FERM,  ESP  and 
NEC  codes  have  aU  been  applied  by  the  Ad¬ 
vanced  Hehcopter  Electromagnetics  (AHE)  con¬ 
sortium  for  analysis  of  antennas  on  helicopters. 
In  this  paper,  helicopter  antennas  are  analyzed 
using  the  FERM  code. 

II.  CODE  MODIFICATIONS 

The  FERM  code  is  written  in  a  modular  struc¬ 
ture,  allowing  it  to  be  highly  portable.  The 
more  computational  intensive  portions  of  the 
code  can  therefore  be  transported  to  comput¬ 
ers  with  larger  memory.  The  modular  form  of 
the  FERM  code  is  a  distinct  advantage.  The 
code  was  originally  installed  on  a  VAX  3100 
workstation.  The  VAX  workstation  required 
too  much  time  to  simulate  antenna  radiation 
problems  related  to  hehcopter  apphcations.  In 
fact,  if  the  number  of  unknowns  reached  1800, 
the  VAX  was  incapable  of  completing  the  com¬ 
putation  due  to  memory  hmitations.  The  VAX 
did  not  have  the  memory  needed  for  the  seg¬ 
mentation  most  geometries  required;  therefore, 
the  code  was  transported  to  an  IBM  RISC/6000 
350  workstation. 

At  first,  the  more  computationally  intensive 
programs  were  transported  to  the  IBM  RISC 
workstation.  Since  antenna  radiation  studies 
were  the  main  areas  of  interest,  the  programs 
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EFIE2,  EFIE3,  and  EFIE4V  were  the  first  to 
be  transferred.  The  program  EFIE2  numeric 
cally  evaluates  the  EFIE’s  for  the  desired  ge¬ 
ometry  and  creates  the  impedance  or  Z-matrix. 
The  program  EFIE3  implements  T  {/decomposi¬ 
tion  to  decompose  the  Z-matrix.  EFIE4V  solves 
for  the  current  distribution  using  the  decom¬ 
posed  Z-matrix  and  the  voltage  sources  stored 
in  the  primary  data  file.  The  execution  of 
EFIEl,  EFIE5V,  and  EFIE5R  remained  on  the 
VAX  workstation.  The  program  EFIEl  con¬ 
verts  the  user  geometry  directives  into  an  inter¬ 
nal  format  for  the  numerical  processing  accom¬ 
plished  via  EFIE2,  EFIE3,  and  EFIE4V.  Even¬ 
tually,  EFIEl  was  migrated  to  the  IBM  RISC 
to  avoid  the  problems  associated  with  transfer¬ 
ring  data  from  VMS  system  to  a  UNIX  system. 
EFIE5V  and  EFIE5R  interface  with  the  graph¬ 
ics  package  DISSPLA.  Since  DISSPLA  is  li¬ 
censed  strictly  to  the  VAX  workstation,  these 
two  routines  must  be  executed  on  the  VAX. 
These  last  two  programs  output  the  antenna  ra¬ 
diation  patterns. 

The  routines  EFIE2,  EFIE3,  and  EFIE4V 
required  minor  modifications  for  execution  on 
the  IBM  RISC.  They  consisted  of  disengaging 
the  subroutine  used  to  calculate  the  CPU  time, 
modifying  the  OPEN  directives  for  temporary 
data  files,  and  altering  the  format  of  the  data 
files  created  by  the  FERM  code.  The  IBM  RISC 
is  capable  of  listing  the  real-time  and  CPU  time 
used  during  the  execution  of  a  program.  There¬ 
fore,  the  subroutine  CLOCKO,  used  to  deter¬ 
mine  the  execution  time  of  a  program  on  a  VAX 
workstation,  was  eliminated. 

AU  calls  to  temporary  or  scratch  files  were 
modified.  VAX/ VMS  Fortran  requires  that 
when  a  SCRATCH  file  is  created  or  opened  the 
initial  size  of  the  file  must  be  specified  using 
the  INITIALSIZE  specifier  [4].  Otherwise,  no 
memory  allocation  is  made.  VS  Fortran  does 
not  recognize  the  specifier  INITIALSIZE  and  it 
was  deleted  from  the  OPEN  statement  for  use 
on  the  IBM  RISC  [5]. 

The  internally-formatted  data  files  created  by 
the  FERM  are  unformatted  and  are  not  read¬ 
able  text  files.  These  files  cannot  traverse  the 
platform  from  a  VAX  workstation  to  an  IBM 


workstation.  Data  was  lost  during  the  use  of 
the  File  Transfer  Protocol  (FTP)  to  transport 
the  data  from  EFIE4V  (IBM  workstation)  to 
EFIE5R  (VAX  workstation).  To  circumvent 
this  problem,  the  data  files  written  to  and  read 
from  (STORAGE,  RESMAT,  and  RESULT)  by 
the  subprograms  had  to  be  formatted  and  re¬ 
quired  commas  between  each  piece  of  data  to 
ensure  no  data  was  lost.  This  involved  rewrit¬ 
ing  aU  the  format  statements  in  EFIEl,  EFIE2, 
EFIE3,  and  EFIE4V,  as  weU  as  the  subroutines 
used  in  these  programs  such  as  MISC  and  EFIE. 

DRAW,  a  graphics  package  developed  at  Ari¬ 
zona  State  University  to  prepare  rectangular 
and  polar  plots  of  radiation  patterns,  has  re¬ 
placed  DISSPLA.  DISSPLA  is  licensed  strictly 
for  execution  on  the  VAX  3100  workstation. 
DRAW  is  executed  on  a  UNIX  workstation.  In 
addition  to  the  ease  of  interface  with  the  FERM 
code  on  the  IBM  RISC,  DRAW  produces  higher 
quality  rectangular  and  polar  plots  than  DIS¬ 
SPLA.  DISSPLA  was  also  used  for  the  ini¬ 
tial  viewing  of  aircraft  geometry.  On  a  UNIX 
based  workstation  this  is  now  done  using  GE- 
OMVIEW  [9],  which  is  available  via  an  anony¬ 
mous  ftp  site. 

In  addition  to  eliminating  DISSPLA  from  use 
as  a  graphics  package,  the  Pattern  Analysis  and 
Database  (PAD)  was  dismantled.  PAD  is  not 
user-friendly.  It  is  capable  of  organizing  and 
manipulating  data  and  storing  numerous  radi¬ 
ation  patterns;  however,  for  antenna  radiation 
prediction,  aU  that  is  needed  is  the  computation 
and  storage  of  radiation  patterns.  PAD  was  re¬ 
placed  by  PLOTR,  a  subroutine  created  specifi¬ 
cally  to  process  radiation  pattern  computations. 
PLOTR  uses  the  SPATTERN  output  data  file 
created  by  EFIE5R  and  calculates  normalized 
or  absolute  antenna  gain. 

III.  GEOM  INTERFACE 

Super  3D  [6]  is  a  popular  interactive  graphics 
package  capable  of  constructing  solid  surface  he- 
hcopter  models.  This  can  be  done  by  employ¬ 
ing  either  a  mouse  or  an  input  data  file.  The 
geometry  can  be  rotated  or  scaled,  thus  allow- 
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ing  modifications  to  be  accomplished  easily  and 
accurately.  Such  a  graphics  package  is  impor¬ 
tant  in  the  analysis  and  design  of  antennas  on 
helicopters.  Super  3D  is  used  for  constructing 
the  helicopter  geometry  and  attaching  the  var¬ 
ious  antennas  to  the  fuselage.  The  solid  sur¬ 
face  geometry  is  then  exported  in  ASCII  for¬ 
mat.  GEOM,  a  geometry  interface  between  the 
Super  3D  and  the  electromagnetics  codes,  trans¬ 
lates  the  solid  surface  model  to  either  a  wire  grid 
model  for  NEC  or  surface  patch  model  for  ESP 

[7],  [8]. 

The  GEOM  program  was  originally  con¬ 
structed  to  allow  interaction  between  the  NEC 
and  ESP  codes  and  the  Super  3D.  However, 
with  the  increasing  use  of  the  FERM  code,  by 
the  AHE  consortium,  for  antenna  radiation  pat¬ 
tern  predictions  and  input  impedance  studies  of 
HE  antennas,  a  link  between  the  FERM  and 
Super  3D  was  also  necessary.  A  subroutine,  to 
convert  geometry  files  from  the  Super  3D  format 
to  the  FERM  format,  was  also  incorporated  into 
the  existing  GEOM  program. 

The  ESP  code  is  a  surface  patch  code,  similar 
to  the  FERM.  However,  the  method  it  employs 
in  modeling  geometries  is  much  different.  The 
ESP  creates  a  surface  by  piecing  together  three- 
dimensional  patches.  Meanwhile,  the  FERM 
uses  a  series  of  lines  to  create  a  two-dimensional 
polygon  representing  a  cross-section  of  the  ge¬ 
ometry.  This  two-dimensional  polygon  is  then 
expanded  along  the  third  dimension  to  join  with 
another  cross-sectional  polygon.  The  circum¬ 
ference  of  the  polygon  may  also  be  scaled  as  it 
is  stretched  to  form  conical  shaped  structures. 
The  FERM  generates  the  triangular  patches  via 
the  combination  of  the  segmentation  used  be¬ 
tween  two  polygons  and  the  segmentation  used 
for  the  lines  of  the  polygon. 

The  FERM’s  method  of  constructing  geome¬ 
tries  is  very  efficient  for  rectangular,  cyhndri- 
cal,  and  spherical  structures,  but  diminishes  as 
the  geometries  become  more  irregular  in  shape. 
This  allows  the  FERM  to  be  highly  suited  for 
modeling  aircraft  fuselages,  but  lacking  in  the 
construction  of  the  stabilizers  and  wings. 

This  modeling  technique  also  eases  the  bur¬ 
dens  encountered  during  segmentation.  As 


mentioned  before,  when  the  FERM  segments  a 
surface,  it  relies  upon  the  number  of  segments 
per  line  comprising  the  cross-sectional  polygon 
and  the  number  of  sections  between  the  poly¬ 
gons.  Therefore,  the  number  of  patches  in  one 
local  area  can  be  increased  without  affecting  the 
segmentation  at  other  areas.  For  example,  the 
number  of  patches  along  the  section  of  the  fuse¬ 
lage  where  the  antenna  is  attached  may  be  in¬ 
creased  while  the  segmentation  along  nose  and 
tail  of  the  aircraft  remains  unaffected. 

During  the  conversion  of  a  Super  3D  geom¬ 
etry  file  to  a  FERM  format,  sufficient  cross- 
sections  of  the  geometry  must  be  taken  to  en¬ 
sure  that  every  contour  of  the  surface  is  closely 
followed.  The  Super  3D  cross-section  coordi¬ 
nates  are  read  into  the  FERM  and  used  to  con¬ 
struct  the  lines  of  the  polygon.  The  number  of 
segments  per  line  is  input  by  the  user.  After 
a  cross-sectional  polygon  is  built,  it  is  joined  to 
the  previous  one  to  build  the  fuselage.  However, 
care  must  be  taken  when  joining  two  polygons. 
The  number  of  segments  per  line  in  a  polygon 
must  be  identical  for  two  polygons  to  be  at¬ 
tached.  Because  of  this,  it  was  decided  that 
the  segmentation  of  each  line,  supplied  by  the 
user,  remain  constant  throughout  the  FERM 
data  file.  The  number  of  segments  per  line  of 
a  cross-sectional  polygon  determines  the  total 
segmentation  around  the  circumference  of  the 
fuselage.  Likewise,  the  total  number  of  sections 
used  between  each  polygon  determines  the  to¬ 
tal  segmentation  of  the  fuselage  length.  If  an 
increase  in  the  segmentation  of  one  area  is  de¬ 
sired,  the  number  of  segments  per  line  in  that 
local  area  may  be  easily  modified  by  the  user  to 
suit  the  application. 

The  stabilizers  and  wings  of  an  aircraft  are 
easily  modeled  with  the  ESP;  they  are  simply 
two-dimensional  plates.  However,  the  FERM  is 
mostly  suited  toward  three-dimensional  applica¬ 
tions.  The  wings  and  stabilizers  are  much  more 
difficult  to  model  than  the  fuselage.  The  FERM 
generates  plates  in  three  ways.  If  the  PLATE 
directive  is  employed,  only  two  lines  with  a  com¬ 
mon  end  point  are  needed.  The  FERM  then 
creates  a  plate  by  generating  two  lines  of  the 
same  length  as  the  original  two  lines  and  plac- 
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ing  them  parallel.  Thus,  the  PLATE  command 
is  limited  to  generating  squares  or  rhomboids. 
If  a  plate  has  more  than  four  sides,  or  has  two 
opposite  sides  that  are  not  parallel,  then  this 
directive  cannot  be  used. 

Implementing  the  TRIA  directive  rather  than 
the  PLATE  directive  provides  a  more  flexible 
method.  The  TRIA  (short  for  triangle)  com¬ 
mand  also  requires  two  lines  with  a  common 
end  point,  but  generates  a  triangle  rather  than  a 
rhombus.  Moreover,  any  two-dimensional  plate 
can  be  constructed  from  a  coUage  of  triangular 
patches.  Therefore,  this  method  provides  the 
most  flexibility  in  generating  arbitrary  shapes 
with  minimal  difiiculty. 


(a) 


(b) 


Figure  1:  The  use  of  the  PRISM  directive  to 
produce  (a)  a  rectangular  plate  and  (b)  a  ta¬ 
pered  wing. 

Another  technique  employed  by  the  PERM 
to  build  two-dimensional  plates  follows  along 
the  same  lines  as  the  construction  of  the  three- 
dimensional  surfaces.  However,  it  is  limited  to 


working  on  plates  with  four  sides  where  at  least 
one  pair  of  opposite  sides  are  parallel.  The 
PRISM  directive  takes  a  polygon  and  deposits 
copies  of  it  at  specified  intervals  along  the  di¬ 
rection  toward  the  desired  ending  point.  A  line 
may  be  redefined  as  a  polygon  and  used  to  gen¬ 
erate  a  rectangular  plate.  In  addition,  the  line’s 
length  may  be  scaled  as  it  is  moved  from  in¬ 
terval  to  interval  so  as  to  make  tapered  plates. 
Figure  1(a)  exhibits  the  use  of  the  PRISM  com¬ 
mand  for  generating  a  rectangular  plate.  The 
original  line,  represented  by  the  solid  line,  is 
separated  into  two  segments  to  provide  a  width 
of  two  segments.  The  dashed  lines  represent  the 
copies  of  the  original  line  placed  so  that  they 
divide  the  length  of  the  plate  into  six  sections. 
While  the  segmentation  along  the  width  of  the 
plate  is  determined  by  the  segmentation  of  the 
original  line,  the  segmentation  along  the  length 
is  accomplished  via  a  numerical  argument  in  the 
PRISM  command.  The  dotted  lines  represent 
the  two  other  sides  of  the  plate  created  by  im¬ 
plementing  PRISM.  In  Figure  1(b)  a  model  of 
a  wing  is  depicted.  This  structure  was  also  cre¬ 
ated  using  the  PRISM  directive,  but  the  length 
of  the  original  line  was  scaled  by  a  value  less 
than  unity  to  produce  a  narrowing  width  at  each 
section  and  yield  a  tapered  shape. 

The  Super  3D  format  does  not  diff'erentiate 
between  patches  used  for  wings  or  patches  used 
for  the  fuselage.  However,  the  GEOM  subrou¬ 
tine  must  decide  whether  the  set  of  data  it  is 
given  represents  a  cross-sectional  polygon  or  a 
two-dimensional  plate.  If  it  is  a  cross-sectional 
polygon,  the  lines  must  be  segmented,  the  poly¬ 
gon  must  be  built  and  joined  to  the  previous 
polygon,  if  one  exists.  If  a  two-dimensional 
plate  must  be  constructed,  then  the  lines  build¬ 
ing  the  plate  must  be  segmented  and  the  plate 
divided  into  as  many  triangles  as  needed  to 
cover  the  surface  area.  To  diff'erentiate  between 
a  cross-sectional  polygon  and  a  two-dimensional 
plate,  the  Super  3D  data  file  must  be  modi¬ 
fied.  The  lines  of  data  representing  the  two- 
dimensional  plate  must  use  a  “plate”  directive 
rather  than  the  “polygon”  directive  normally 
used.  This  alteration  is  a  necessity  for  the 
GEOM  subroutine  to  perform  correctly.  If  a 
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line  is  encountered,  separate  from  other  lines  or 
polygons,  then  it  is  modeled  as  a  wire  monopole 
antenna. 

The  completion  of  the  GEOM  subroutine  to 
transfer  Super  3D  data  files  to  FERM  data  file 
marks  a  milestone.  Now  the  identical  geometry 
analyzed  by  the  NEC  and  ESP  can  be  studied 
with  the  FERM  without  inaccuracies  resulting 
from  the  reconstruction  of  the  structure. 

IV.  RESULTS 

A.  PATTERN  PREDICTION 

Figures  2  and  3  illustrate  the  geometry  and 
overall  dimensions  of  the  Blackhawk  helicopter, 
and  the  FERM  model  used  in  this  analysis, 
respectively.  A  loop  or  towel-bar  antenna  is 
mounted  beneath  the  tail  of  the  Blackhawk  he¬ 
licopter.  The  loop  had  a  length  of  3.97  meters 
and  a  height  of  0.4  meters.  Because  a  strip 
antenna  was  implemented,  a  width  of  4.0  cm, 
equivalent  to  a  wire  radius  of  1.0  cm,  was  used 
[10].  The  loop  consists  of  only  three  strips  with 
the  metallic  body  of  the  helicopter  comprising 
the  fourth  side.  The  antenna  is  fed  along  the 
strip  closest  to  the  front  of  the  helicopter. 

There  are  two  problems  associated  with  the 
FERM  code  in  modeling  wire  loop  antennas 
such  as  the  one  illustrated  in  Figure  3.  The  first 
is  that  two  wires  must  be  coUinear  or  connected 
end-to-end.  In  generating  a  loop  antenna,  the 
wires  cannot  be  placed  end-to-end  at  the  90® 
bends.  The  second  problem  is  that  wires  can¬ 
not  be  electrically  attached  to  surfaces  in  the 
current  version  of  the  FERM  code.  To  solve  the 
first  problem,  the  loop  is  created  using  a  strip 
equivalent  model.  The  strip  can  be  electrically 
connected  at  right  angles  while  the  wires  can¬ 
not.  Figure  4  exhibits  the  differences  between 
the  wire  and  strip  models  at  the  corner  of  the 
loop  antenna.  The  strip  can  also  be  electrically 
attached  to  the  helicopter  fuselage  if  there  are 
two  common  nodes  between  the  strip  mesh  and 
the  fuselage  mesh.  This  was  not  possible  to  do 
for  the  loop  antenna  which  was  modeled  as  a 
strip  due  to  the  changing  size  of  the  helicopter 
tail,  and  so  the  edges  of  the  strip  are  just  placed 


near  the  surface  of  the  helicopter  fuselage.  This 
is  not  expected  to  influence  much  the  computa¬ 
tion  of  the  antenna  radiation  patterns. 

Figure  5  displays  the  radiation  patterns  of 
the  loop  antenna  mounted  on  the  Blackhawk 
helicopter  at  10  MHz.  FERM  and  NEC  com¬ 
puted  radiation  patterns  for  the  co-polarized 
yaw  plane  are  compared.  The  patterns  obtained 
by  the  two  codes  were  normalized  to  0  dB  based 
on  their  individual  maxima.  The  FERM  ex¬ 
hibits  a  magnitude  approximately  2  dB  greater 
than  that  of  the  NEC  around  the  ^  =  90®  and 
<j)  =  270®.  However,  the  results  from  the  two 
codes  agree  weU  in  the  remaining  portions  of 
the  pattern.  In  Figure  6,  the  co-polarized  pitch 
plane  radiation  patterns  are  shown.  Unlike  the 
yaw  plane,  the  agreement  between  the  NEC  and 
FERM  in  the  pitch  plane  is  not  as  good.  The 
co-polarized  roU  plane  is  shown  in  Figure  7.  The 
best  agreement  between  the  two  codes  is  seen  in 
this  plane. 


Side  View  Top  View 


2J«2m 


Figure  2:  Geometry  of  the  Blackhawk  heli¬ 
copter. 

B.  INPUT  IMPEDANCE  PREDIC¬ 
TION 

For  a  strip  antenna  to  be  electrically  connected 
to  the  helicopter  fuselage  or  any  other  geometry, 
it  must  have  a  common  edge  with  that  geome¬ 
try;  i.e.,  a  multiple  edge  between  the  strip  an¬ 
tenna  and  the  ground  plane  (or  helicopter  fuse¬ 
lage)  must  exist  for  the  two  to  be  electrically 
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Figure  3:  Geometry  of  the  4.0  m  towel-bar  an¬ 
tenna  mounted  on  the  segmented  Blackhawk  he¬ 
licopter. 


(a)  (b) 

Figure  4:  The  modeling  of  a  right  angle  of  a 
loop  antenna  using  (a)  wires  and  (b)  strips. 


Figure  5:  The  co-polarized  yaw  plane  radiation 
pattern  of  the  loop  antenna  on  the  Blackhawk 
helicopter  model  at  10  MHz.  There  is  no  rotor 
blade  structure. 

0 


Bottom 


Figure  6:  The  co-polarized  pitch  plane  radiation 
pattern  of  the  loop  antenna  on  the  Blackhawk 
helicopter  model  at  10  MHz.  There  is  no  rotor 
blade  structure. 
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Figure  7:  The  co-polarized  roll  plaue  radiation 
pattern  of  the  loop  antenna  on  the  Blackhawk 
helicopter  model  at  10  MHz.  There  is  no  rotor 
blade  structure. 


connected. 

To  illustrate  how  a  strip  antenna  is  prop¬ 
erly  connected  to  a  ground,  a  simpler  geometry 
is  considered.  The  geometry  is  that  of  a  A/4 
monopole  on  a  1  x  lA  square  ground  plane  at 
100  MHz.  Of  interest  in  this  problem  is  the 
calculation  of  the  input  impedance  of  the  an¬ 
tenna.  The  geometry  of  the  strip  antenna  and 
the  way  it  is  properly  attached  to  the  ground 
plate  is  illustrated  in  Figure  8.  The  thickness 
of  the  strip  was  chosen  four  times  the  wire  ra¬ 
dius  (w=20  mm).  To  connect  the  strip  to  the 
ground  plate  the  two  must  have  a  common  edge 
along  a  triangular  patch.  The  strip  in  this  case 
was  oriented  in  the  xz-plane,  and  therefore,  the 
discretization  of  the  plate  in  the  x  direction  was 
made  nonuniform.  This  was  required  because 
the  grid  size  in  the  x-direction  near  the  middle 
of  the  plate  (where  the  strip  is  attached)  must 
be  equcd  to  the  width  of  the  strip.  This  is  illus¬ 
trated  in  Figure  8. 

Using  the  simple  geometry  of  the  monopole 
on  a  squaxe  ground  plate  the  input  impedance 
of  the  antenna  was  calculated  in  the  frequency 
range  from  10  to  100  MHz.  The  FERM  com¬ 
puted  input  resistance  is  compared  with  similar 
predictions  using  the  NEC  code  in  Figure  9,  and 
the  antenna  input  reactance  in  Figure  10.  The 
results  of  having  the  strip  antenna  on  top  of  the 


Figure  8:  An  example  of  a  strip  antenna  con¬ 
nected  to  a  ground  plate  for  FERM  modeling. 


plate  without  being  electrically  connected  to  the 
plate  are  also  included.  As  illustrated  by  both 
figures,  there  is  only  a  few  Ohms  difference  be¬ 
tween  the  FERM  predictions  with  the  properly 
connected  antenna  and  the  NEC  results. 

It  is  expected  that  if  the  monopole  antenna 
is  properly  attached  to  the  helicopter  fuselage 
then  similar  input  impedance  results  will  be  ob¬ 
tained,  especially  at  HF  frequencies.  It  was  pre¬ 
viously  demonstrated  that  the  current  distribu¬ 
tion  at  HF  frequencies  is  localized  in  the  vicinity 
of  the  antenna,  and  that  the  platform  used  to 
mount  the  antenna  is  not  expected  to  signifi¬ 
cantly  affect  the  input  impedance  calculations 
[11]. 

C.  EXECUTION  TIMES  OF  THE 
FERM  ON  DIFFERENT  COM¬ 
PUTERS 

A  comparison  of  the  execution  times  of  the 
FERM  code  subroutines  was  accomplished  to 
determine  whether  the  IBM  RISC  6000  is 
more  efficient  than  the  VAX  3100.  Since  the 
VAX  3100  workstation  cannot  handle  the  large 
scratch  files  generated  for  geometries  with  more 
than  1500  unknowns,  the  scope  of  the  compari¬ 
son  was  bounded  by  this  upper  limit.  The  four 
main  numerical  processing  subroutines  used  in 
predicting  antenna  radiation  patterns  were  con¬ 
sidered:  EFIEl,  EFIE2,  EFIE3,  and  EFIE4V. 
The  results  of  this  test  are  displayed  in  Table  1. 


Table  1:  Execution  times  for  PERM  subroutines 
on  various  computers. 


Figure  9:  NEC  and  FERM  computed  input  re¬ 
sistance  of  a  strip  antenna  on  a  finite-size  con¬ 
ducting  ground  plate,  in  the  frequency  range 
from  10  to  100  MHz. 
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Figure  10:  NEC  and  FERM  computed  input  re¬ 
actance  of  a  strip  antenna  on  a  finite-size  con¬ 
ducting  ground  plate,  in  the  frequency  range 
from  10  to  100  MHz. 


Comp. 

Type 

No. 

of 

Unkn. 

Execution  Time 
(seconds) 

EFIEl 

EFIE2 

EFIE3 

EFIE4V 

VAX 

299 

1.0 

209.0 

28.0 

11.0 

IBM 

299 

0.05 

0.72 

0.440 

3.26 

VAX 

517 

3.0 

840.0 

103.0 

96.0 

IBM 

517 

1.3 

150.2 

81.6 

73.8 

VAX 

1201 

10.0 

4570 

842 

375 

IBM 

1201 

4.7 

1661.7 

818.1 

300.2 

For  the  first  two  entries  in  the  table,  a  A/4 
wire  monopole  located  over  a  perfectly  con¬ 
ducting  ground  plane  was  used  for  comparison. 
At  the  test  frequency  of  50  MHz,  the  physi¬ 
cal  length  of  the  antenna  was  1.5  meters  while 
the  sides  of  the  ground  plane  were  6.0  meters 
in  length.  The  diameter  of  the  antenna  was 
1.5  cm.  The  gap  between  the  antenna  and  the 
ground  plane  was  matched  to  the  diameter  of 
the  wire.  The  plate  is  segmented  into  100  tri¬ 
angular  pairs  (ten  segments  per  side).  The  wire 
antenna  is  divided  into  ten  segments  with  the 
voltage  excitation  placed  between  the  first  and 
second  antenna  segments  closest  to  the  plate. 
This  configuration  spawned  299  unknowns. 

There  was  little  difference  between  execu¬ 
tion  times  on  the  VAX  and  the  IBM  RISC  for 
the  first  (EFIEl),  third  (EFIE3),  and  fourth 
(EFIE4V)  subroutines.  The  VAX  workstation 
produced  results  within  30  seconds  of  the  IBM 
RISC  for  these  subroutines.  There  was  a  sub¬ 
stantial  difference  in  run  times  for  the  second 
subroutine,  EFIE2.  This  routine  evaluates  the 
electric  field  integral  equation  for  the  input 
geometry,  creates  an  impedance  matrix,  and 
writes  the  matrix  to  a  data  file.  This  com¬ 
putation  is  of  order  N^,  where  N  is  the  num¬ 
ber  of  unknowns.  The  IBM  RISC  performed 
this  task  in  under  1  second  of  CPU  time,  while 
the  VAX  required  approximately  200  seconds 
of  CPU  time.  Since  the  number  of  unknowns  is 
relatively  small,  the  amount  of  time  needed  to 


write  to  and  read  from  the  scratch  file  TEMP 
is  not  significant.  Therefore,  the  difference  in 
execution  times  for  the  second  subroutine  can 
be  attributed  to  the  IBM  RISC’s  more  efficient 
numerical  evaluation  of  the  electric  field  integral 
equations. 

A  block  helicopter  modeled  with  a  strip  an¬ 
tenna,  shown  in  Figure  11(b),  provided  the  next 
set  of  unknowns  for  comparison.  The  test  fre¬ 
quency  for  this  geometry  was  established  as  70 
MHz.  The  block  helicopter  length  (4.0  m)  was 
partitioned  into  10  segments.  The  A/4  (1.07  m) 
antenna  also  had  10  segments  with  a  1.0  volt 
excitation  located  between  the  first  and  second 
segments  closest  to  the  fuselage.  This  geometry 
yielded  517  unknowns  for  processing  on  the  two 
workstations. 

This  portion  of  the  examination  was  very 
similar  to  the  test  runs  with  299  unknowns. 
The  VAX  workstation  had  execution  times  22- 
23  seconds  slower  than  those  of  the  IBM  RISC 
for  the  EFIE3  and  EFIE4V  subroutines.  EFIE2 
took  5.6  times  as  long  to  execute  on  the  VAX 
as  it  did  on  the  IBM  RISC.  This  case  differed 
from  the  first  in  that  EFIE2  had  to  use  disk 
memory  to  complete  the  numerical  evaluation 
on  the  VAX.  This  process  slows  the  execution 
time  substantially  on  the  VAX.  More  on  this 
subject  will  be  seen  as  the  number  of  unknowns 
increases  towards  the  upper  limit. 

The  final  geometry  used  for  this  comparison 
was  a  modified  version  of  the  Blackhawk  heli¬ 
copter  in  Figure  3.  The  segmentation  along  the 
length  of  the  helicopter  was  decreased  as  was 
the  segmentation  in  the  wings  and  stabilizers. 
This  was  done  so  that  the  number  of  unknowns 
would  remain  under  1500.  This  modified  heli¬ 
copter  geometry  produced  1201  unknowns. 

As  the  number  of  unknowns  was  quadru¬ 
pled  from  299  to  1201,  the  difference  in  CPU 
times  for  the  EFIE3  subroutine  remained  con¬ 
stant.  EFIE3  is  the  subroutine  that  implements 
L  U  decomposition  to  decompose  the  impedance 
matrix  created  by  the  EFIE2  subroutine.  Sim¬ 
ilarly,  the  differences  in  execution  times  for  the 
EFIEl  and  EFIE4V  increased  very  slightly  with 
the  increase  in  the  number  of  unknowns.  The 
5.3  second  difference  between  the  VAX’s  run 


time  and  the  IBM  RISC’s  run  time  for  the 
EFIEl  subroutine  is  insignificant.  As  the  num¬ 
ber  of  unknowns  was  increased  from  517  to 
1201,  the  execution  time  on  the  IBM  RISC 
for  the  second  subroutine  increased  by  more 
than  an  order  of  magnitude.  Meanwhile,  the 
VAX’s  execution  time  for  the  same  subroutine 
only  increased  by  a  factor  of  five.  Clearly,  for 
the  case  of  1201  unknowns,  the  difference  in 
CPU  time  between  VAX  and  IBM  RISC  in¬ 
creased,  but  their  ratio  is  decreasing.  This 
suggests  that  as  the  number  of  unknowns  in¬ 
creases,  the  efficiency  of  the  VAX  workstation 
approaches  that  of  the  IBM  RISC  workstation. 
However,  this  conclusion  cannot  be  completely 
verified  because  of  the  memory  limitations  on 
the  VAX  3100  workstation.  For  geometries  with 
less  than  1200  unknowns,  the  IBM  RISC  6000  is 
a  better  computer  platform  for  execution  of  the 
FERM  code  because  it  operates  more  efficiently. 
Yet,  for  geometries  with  a  larger  number  of 
unknowns,  the  IBM  RISC  6000  is  stiU  better, 
because  of  the  substantially  larger  amount  of 
memory  required  by  the  temporary  scratch  files 
written  during  the  numerical  evaluation  of  the 
electric  field  integral  equations  in  the  subroutine 
EFIE2. 

V.  CONCLUSIONS 

The  Finite  Element  Radiation  Model  code  is 
a  Method  of  Moment  solver  used  primarily  for 
radar  cross  section  prediction  and  antenna  mod¬ 
eling.  The  code  has  been  applied  in  this  paper 
to  predict  the  antenna  radiation  patterns  of  a 
loop  mounted  on  the  Blackhawk  helicopter  and 
the  input  impedance  of  a  monopole  antenna  on 
a  ground  plate.  The  results  have  shown  that 
when  modeling  a  loop  or  towel-bar  antenna  with 
the  FERM  code,  the  antenna  must  be  com¬ 
prised  of  thin  strips  rather  than  wires.  This  en¬ 
sures  that  the  corners  of  the  loop  are  connected 
both  physically  and  electrically.  The  width  of 
the  strips  was  taken  to  be  four  times  the  equiv¬ 
alent  radius  of  a  wire  with  circular  cross-section 
[10].  However,  there  are  some  difficulties  in 
properly  attaching  the  loop  edges  to  the  heli- 
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Figure  11:  Antenna  geometries  on  the  block  he¬ 
licopter  model:  (a)  wire  antenna,  (b)  strip  an¬ 
tenna  with  the  width  along  the  x-axis,  and  (c) 
strip  antenna  with  the  width  along  the  y-axis. 

copter  surface  due  to  the  complexity  of  the  fi¬ 
nite  element  mesh.  To  alleviate  this  problem 
the  loop  antenna  was  placed  near  the  helicopter 
surface.  The  radiation  patterns  resulting  from 
this  analysis  were  compared  with  similar  ones 
obtained  using  the  NEC  code.  For  antenna  in¬ 
put  impedance  calculations  a  simpler  structure 
was  considered,  that  of  a  monopole  on  a  ground 
plate.  It  is  demonstrated  that  accurate  input 
impedance  results  can  be  obtained  when  the  an¬ 
tenna  is  properly  grounded  on  the  plate. 
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Abstract 

The  Finite  Element  Method  with  second  order 
Absorbing  Boundary  Condition  is  a  recently 
developed  computational  technique  that  finds 
use  in  antenna  design  and  Electromagnetic 
Compatibility  simulation.  To  determine  the 
accuracy  of  this  new  procedure,  the  problem  of 
aperture  radiation  was  studied.  The  near  zone 
and  aperture  fields  of  a  waveguide  antenna 
computed  using  the  Finite  Element  Method  are 
compared  with  published  data  and  results 
found  using  other  simulation  techniques. 
These  comparisons  show  that  the  new  method 
is  accurate  even  when  the  ABC  boundary 
conforms  to  the  problem  geometry  and  is 

placed  as  near  as  Xflit  to  the  aperture. 

Introduction 

A  central  problem  in  antenna  design  and 
Electromagnetic  Compatibility  (EMC)  is 
computing  the  radiation  from  apertures. 
Apertures  represent  the  key  element  in  many 
radiating  structures  including  open  ended 
waveguides,  flared  horns  and  reflector 
antennas.  In  EMC,  openings  in  system 
cabinets  provide  a  major  source  of  unwanted 
radiation.  This  latter  application  has  become 
increasingly  important  as  the  European 
Economic  Community  and  the  FCC  tighten 
their  EMC  compliance  regulations. 

Aperture  problems  have  been  studied  in 
the  past  [1  and  2].  Typically  the  Method  of 
Moments  (MOM)  is  used  to  solve  for  the  fields 
in  the  aperture  assuming  that  the  aperture  is  cut 
in  an  infinite  ground  plane.  While  this 
provides  an  efficient  method  of  solution,  the 
infinite  ground  plane  assumption  limits  its 
utility.  The  resulting  model  is  appropriate  for 
analyzing  an  array  of  aperture  antennas,  but  not 
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for  computing  the  radiation  from  finite  sized 
enclosures.  Real  world  systems  are  finite  in 
size  and  often  have  irregular  openings  and 
other  objects,  both  metal  and  dielectric,  that 
affect  the  radiation.  For  such  problems  a 
hybrid  approach  can  be  used  [3  and  4]. 
Unfortunately  these  methods  are  only 
appropriate  for  the  specific  geometries 
considered.  To  solve  more  general  problems,  a 
general  purpose  computer  simulation  tool  is 
needed. 

A  technique  that  can  easily  handle 
complex  geometries  and  provides  a  general 
purpose  simulation  procedure  is  the  Finite 
Element  Method  (FEM).  The  FEM  has  been 
used  successfully  in  the  past  to  model  closed 
region  problems.  In  the  FEM,  the  solution 
space  is  broken  into  small  elements  and  the 
field  in  each  element  is  calculated  directly.  The 
resulting  discretization  is  usually  referred  to  as 
the  finite  element  mesh.  For  open  region 
radiation  problems,  the  solution  space  must  be 
bounded.  One  way  to  accomplish  this  is  to 

place  a  377  Cl  impedance  boundary  far  from  the 
radiating  source.  However,  this  approach 
generates  a  very  large  solution  region  and 
requires  correspondingly  large  solution  times. 
Recently  a  new  Absorbing  Boundary  Condition 
(ABC)  has  been  developed  [5  and  10]  that  can 
be  used  to  terminate  the  FEM  mesh  at 
boundaries  very  close  to  the  radiating  structure. 
The  smaller  solution  region  with  this  new 
absorbing  boundary  makes  the  resulting 
computation  more  efficient.  However,  the 
following  question  arises:  How  close  can  the 
new  ABC  be  placed  to  the  radiating  structure 
without  significant  loss  of  accuracy?  To 
answer  this  question,  a  study  was  undertaken 
to  compare  the  results  from  FEM  simulations 
with  measured  data  available  from  the  literature 
and  with  data  found  by  using  other  solution 
techniques. 

While  the  radiation  from  an  opening  in  a 
system  cabinet  is  an  important  problem,  the 


available  data  is  limited.  Therefore,  the 
radiation  from  a  waveguide  antenna  was 
chosen  as  the  test  case  for  this  study.  Using  a 
waveguide  antenna,  it  will  be  shown  that  this 
new  ABC  can  be  made  to  conform  to  the 
problem  geometry  and  can  be  placed  as  close  as 

XHti  from  the  aperture  for  such  geometries. 
Placing  the  ABC  this  close  to  the  radiating 
source  keeps  the  FEM  mesh  small  and  results 
in  an  efficient  solution. 

Accuracy  Statement 

For  the  case  presented  here  where  the  ABC 
is  at  XI2%  the  magnitude  of  the  near  zone 
electric  field  along  the  aperture  center  line  is 
accurate  to  within  ±  0.5  dB  when  compared  to 
measured  data.  Also  the  calculation  of  the 
amplitude  of  the  dominant  mode  in  the  aperture 
of  an  X  band  waveguide  for  the  same  ABC 
spacing  is  accurate  to  within  2  %  in  magnitude 
and  2°  in  phase  when  compared  to  a  method  of 
moments  calculation.  These  numbers  were 
found  to  be  typical  for  the  aperture  radiation 
problems  considered  here.  The  apertures 
considered  here  have  a  maximum  dimension 

that  is  on  the  order  of  1  X,.  As  an  example  of 
the  resources  required  for  a  solution:  for  the 
simulation  of  the  unflanged  waveguide  with  the 

ABC  at  XHn  (data  presented  in  Figure  3)  the 
run  time  for  3  solution  passes  was  less  than 
one  hour  on  an  HP  735  workstation  with  a  99 
MHz.  processor  and  400  MB  memory. 

Overview  of  the  FEM 

To  understand  the  results  presented 
here  it  is  useful  to  first  review  the  basics  of  the 
FEM.  See  Reference  7  for  more  complete 
details  on  the  FEM.  In  the  FEM,  the  solution 
space  is  broken  into  small  regions  called  finite 
elements.  The  elements  used  in 
electromagnetics  are  often  tetrahedral  as 
illustrated  in  Figure  1.  This  element  has  the 
advantage  of  being  both  flexible  for  modeling 
shapes  and  allowing  closed  form  integrations. 
Since  the  solution  space  is  split  into  small 
elements  it  is  very  easy  to  model  geometries 
that  contain  different  materials.  This  is  most 
helpful  when  simulating  printed  circuit 
geometries. 


The  electric  field  in  each  element  is 
approximated  by  using  the  tangential 
components  along  each  edge  as  shown  in 
Figure  1.  This  choice  of  approximation 
functions  eliminates  the  spurious  modes  that 
can  cause  problems  in  finite  element 
simulations  [8].  The  tangential  electric  fields 
are  approximated  by  polynomials  that  contain 
unknown  weighting  coefficients,  Ej.  The 

electric  field  must  satisfy  the  vector  wave 
equation: 

VxVxE  =  E  (1) 

The  vector  wave  equation  leads  to  the 
following  variational  principle: 

F(i)  =  j  (vxif  d  s  -  j  (if  d  s  (2) 

The  polynomial  approximation 
functions  are  substituted  into  (2)  and  the 
resulting  set  of  equations  are  minimized  with 
respect  to  the  unknown  coefficients  by  setting 
the  3F/9Ej  equal  to  0  for  all  i.  While  the 

resulting  matrix  is  generally  large,  it  is  sparse. 
This  matrix  equation  is  solved  to  find  the 
unknowns,  Ep 

In  FEM  simulation,  the  size,  shape  and 
location  of  the  elements  in  the  mesh  is 
important  to  the  accuracy  of  the  solution.  To 
reduce  the  possibility  of  a  poor  mesh  causing 
inaccuracies,  Delaunay  tessellation  and  adaptive 
mesh  refinement  [16]  was  employed.  In 
adaptive  mesh  refinement,  the  solution  is  first 
computed  with  a  coarse  initial  mesh.  The 
approximate  solution  is  then  substituted  into  the 
vector  wave  equation  (1)  to  determine  the  eiTor 
in  each  element.  Next  the  elements  that  contain 
the  largest  errors  are  refined.  The  result  is  that 
the  elements  in  those  parts  of  the  problem 
where  the  errors  are  the  largest  are  reduced  in 
size.  The  entire  procedure  is  repeated  until  the 
desired  accuracy  is  attained.  In  the  CAE  tool 
used  here  the  user  controls  the  definition  of  the 
initial  mesh.  The  user  can  specify  a  seed  value 
for  various  surfaces  in  the  model.  This  seed 
value  is  the  spacing  of  additional  points  added 
to  the  original  problem  description  on  a 
specified  surface.  These  seed  points  become 
the  vertices  of  the  tetrahedral  finite  elements. 


The  Second  Order  ABC 

Absorbing  boundary  conditions 
(ABCs)  are  of  two  general  types.  The  first  is 
the  one  derived  by  Bayliss  and  Turkel  [9]  in 
scalar  form  and  then  by  Peterson  [5]  and  by 
Webb  and  Kanellopoulos  [10]  in  vector  form. 
The  basis  of  the  derivation  of  this  ABC  is  the 
far  field  expansion  or  Wilcox  expansion  [11]. 
The  first  term  in  this  expansion  provides  a  377 

Q.  impedance  boundary  that  must  be  placed 
relatively  far  from  radiating  sources.  The 
second  type  of  ABC  was  derived  by  Engquist 
and  Majda  [12]  in  scalar  form  and  %  Sun  and 
Balanis  [13]  in  vector  form.  This  ABC  is 
derived  by  splitting  the  operator  into  two  first- 
order  operators  and  then  approximating  these 
first-order  operators.  This  allows  waves  to 
travel  in  one  direction  only;  hence,  they  are 
often  called  a  one-way  ABCs. 

In  this  paper  the  second-order  ABC 
presented  in  References  [5  and  10]  is  used.  It 
should  be  noted  that  although  these  ABCs  are 
derived  in  the  reference,  these  papers  are  purely 
theoretical  and  do  not  provide  numerical 
examples  of  its  accuracy.  The  purpose  of  this 
paper  is  to  provide  such  examples. 

This  ABC  is  computationally  efficient 
for  the  following  reasons.  First,  it  includes 
higher  order  terms  compared  to  the  first  order 
ABC.  Second,  the  ABC  may  be  applied  on  a 
conformal  boundary  because  it  absorbs  the 
outgoing  wave  on  a  piecewise  basis.  Finally, 
the  finite  element  shape  functions  themselves 
are  second  order  so  that  the  ABC  boundary 
may  be  placed  very  close  to  the  source  of  the 
radiation. 

The  above  discussion  indicates  that  the  new 
ABC  behaves  very  much  like  the  absorber  used 
in  an  anechoic  chamber.  Therefore,  one  would 
expect  that  the  second  order  ABC  can  be  placed 
very  close  to  the  radiating  structure  and  can  be 
made  to  conform  to  the  geometry.  From 
antenna  theory  it  can  be  shown  that  the  near 
zone  of  an  antenna  can  be  broken  into  two 
separate  regions:  the  reactive  near  zone  and  the 
radiating  near  zone  regions  [14].  The  reactive 
near  zone  is  the  region  closest  to  the  geometry 
where  the  reactive  fields  are  present.  It  is  the 
fields  in  this  region  that  contribute  to  the 
imaginary  component  of  the  input  impedance. 
Since  these  are  reactive  fields  one  does  not 
want  to  place  an  absorber  in  this  region  else  the 


amount  of  energy  stored  is  affected.  On  the 
other  hand,  this  should  not  be  a  concern  in  the 
radiating  near  zone.  Therefore,  one  should  be 
able  to  place  the  ABC  in  the  radiating  near  zone 
of  the  geometry.  Reference  [14]  states  that  the 
boundary  between  the  reactive  and  radiating 

near  zones  may  be  as  small  as  X/2n  from  the 
antenna.  One  would  expect  then  that  the  lower 
limit  on  the  ABC  placement  would  be  on  the 

order  of  X/2k. 

Results  for  an  Unflanged  Waveguide 

The  computed  data  presented  here  was 
solved  using  an  initial  mesh  that  was  seeded 

with  mesh  points  every  Ay  10  over  the  aperture 
opening,  over  the  input  port  and  on  the  ABC 
boundary.  Adaptive  mesh  refinement  was 
applied  twice  to  these  initial  values  to  perform  a 
total  of  three  solution  passes.  This  was 
sufficient  to  ensure  that  the  change  in  the  value 
of  the  reflection  coefficient  from  pass  2  to  pass 
3  was  less  than  0.02. 

A  problem  with  much  of  the  waveguide 
antenna  simulations  presented  in  the  Uterature  is 
that  the  aperture  plane  is  assumed  to  be  an 
infinite  ground  plane  corresponding  to  a 
waveguide  with  an  infinite  flange.  While  this 
assumption  reduces  the  complexity  of  the 
mathematical  model,  it  does  not  reflect  the 
situation  in  the  real  world.  To  determine  the 
optimum  distance  for  placement  of  the  ABC,  it 
is  important  to  simulate  a  geometry  that  does 
not  have  an  infinite  flange.  The  authors  of 
Reference  15  considered  such  a  problem:  they 
studied  the  radiation  from  an  unflanged 
WR3200  rectangular  waveguide.  This  is  the 
type  of  antenna  used  in  large  anechoic 
chambers  for  susceptibility  testing.  For  the  case 
considered  in  Reference  15  the  antenna  has  a 
net  input  power  of  44.8  watts  at  a  frequency  of 
275  MHz.  The  FEM  simulation  model  that  was 
used  is  shown  in  Figure  2.  This  is  a  simple 
unflanged  waveguide  with  an  input  port  at  one 
end.  The  waveguide  is  enclosed  in  a 
rectangular  box  on  which  is  placed  the  ABC. 
The  spacing  D  was  varied  and  the  data 
compared  with  the  measured  values.  The 
results  are  shown  in  Figure  3.  The  data 
presented  is  the  magnitude  of  the  electric  field 
along  the  x=y=0  line  for  a  variety  of  values  of 
z  in  the  antenna's  near  zone.  At  the  frequency 


of  operation  the  field  points  for  the  smaller 

spacings  are  less  than  1  A,  from  the  opening.  It 
should  be  noted  that  these  field  points  are 
outside  of  the  FEM  mesh  for  the  cases 
presented  here.  The  values  of  the  fields 
determined  on  the  boundary  are  used  to 
calculate  the  fields  at  points  outside  of  the  mesh 
(see  [6]  for  details).  From  the  data  shown  in 

Figure  3,  it  can  be  seen  that  for  D  =  A,/4  and 

Xy27C  the  differences  between  measured  and 
computed  results  are  less  than  +  0.5  dB.  The 

X./10  case,  though,  is  not  as  accurate.  These 
results  are  an  additional  0.5  dB  smaller  than  the 

A/2rt  data  when  compared  to  the  measured.  In 
fact  a  case  could  be  made  that  at  this  spacing 
the  results  are  still  very  good.  The  FEM 
meshes  for  the  3  spacings  had  approximately 

the  same  number  of  unknowns  (within  *10% 
of  each  other).  The  main  difference  between 
these  three  simulations  was  the  size  of  the 
tetrahedra,  not  a  significant  increase  in  the 
number. 

Another  test  of  the  solution  procedure 
was  to  calculate  the  gain  of  the  antenna  using 
the  far  zone  fields  determined  by  the  FEM/ABC 
approach.  The  gain  for  this  antenna  calculated 
using  equation  (2)  in  Reference  15  is  6.9  dB 
and  the  gain  calculated  in  the  simulation  with 

D=X/271  is  6.5  dB.  While  the  gain  equation 
given  in  Reference  15  is  approximate,  it 

nevertheless  indicates  that  a  D=X/2tc  spacing 
can  be  used  to  accurately  model  this  aperture 
problem.  Thus,  the  ABC  can  be  made  to 
conform  to  the  shape  of  the  geometry  and  can 

be  placed  as  close  as  XHn  from  the  aperture 
with  accurate  results  for  these  geometries. 

Using  D=X./4  for  the  spacing  does  not 
significantly  increase  the  accuracy  of  the  results 
but  does  require  a  small  increase  in  mesh  size. 

As  a  final  note,  the  data  for  the  'kUn 

and  X/4  cases  appeared  to  converge  to  accurate 
values  for  the  reflection  coefficient  after  only 
one  adaptive  mesh  refinement,  but  they  were 
run  for  a  second  mesh  refinement  to  be 
consistent  with  the  other  case. 


The  Aperture  Fields  for  Flanged 
Waveguide  Antennas 

As  stated  in  the  introduction,  the 
radiation  from  apertures  has  been  studied  in  the 
past  using  the  method  of  moments  (MOM). 
Thus,  the  results  from  a  MOM  solution  can 
also  be  used  to  verify  the  FEM  simulations. 
Before  discussing  the  specifics  of  this 
approach,  it  should  be  noted  that  this  MOM 
procedure  is  designed  to  analyze  waveguide 
antennas  with  infinite  flanges  only  and  it  is  not 
a  general  purpose  CAE  tool.  In  this 
simulation,  the  aperture  plane  is  assumed  to  be 
an  infinite  ground  plane  and  the  aperture  is 
closed  using  the  equivalence  principle. 
Equivalent  magnetic  currents  are  placed  on 
either  side  of  the  short.  These  magnetic 
currents  are  approximated  by  a  finite  weighted 
sum  of  entire  domain  basis  functions.  The 
unknown  weighting  coefficients  are  determined 
using  a  Galerkin  procedure.  Vector  potential 
modal  functions  were  used  as  basis  functions. 
The  field  approximation  for  the  aperture  shown 
in  Figure  4  is  as  follows  [2]: 

Ea  “  X  [^nm  ®y,nm  y  Dnm  ®x,nm  x]  (3) 

n,m 

where: 


ey.nm  =  sin(S|^)  COs(2^j 

*  ^  ^  (4) 

ex,nm  =  COs(Sf-) 

and  where  80n  is  Neumann's  number.  Since 
the  MOM  simulation  assumes  an  infinite 
flange,  a  small  ground  plane  was  placed  around 
the  aperture  in  the  FEM  models.  The  ground 

plane  extends  A,/4  from  each  edge  of  the 
aperture  in  the  ±x  and  dty  directions.  The  ABC 

is  rectangular  in  shape.  It  is  located  at  D  =  X/27C 
from  the  aperture  and  from  the  edges  of  the 
ground.  The  apertures  and  the  ABC  were 
seeded  in  the  same  way  as  before  and  a  total  of 
2  passes  were  run  for  all  of  the  simulations  in 
this  section. 

Due  to  the  orthogonality  of  the  basis 
functions  of  (4),  it  is  a  simple  task  to  calculate 
the  coefficients  Cj^jj^  and  from  the  FEM 
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data.  To  compute  or  the  dot 

product  of  the  fields  in  the  aperture  are  taken 
with  Cy  for  C^ni  or  with  ex,nm^  for 
Djjjjj,  and  the  result  is  integrated  numerically 

over  the  aperture.  Consider  the  data  presented 
in  Figures  5  and  6.  The  magnitude  and  phase 
of  the  first  4  modes  in  an  X  band  waveguide 
operated  at  10  GHz  are  shown.  It  was  found 
previously  that  these  four  modes  are  sufficient 
to  accurately  model  the  aperture  field  in  this 
antenna  [2].  As  stated  previously,  the  FEM 
and  MOM  data  for  the  dominant  mode  are 
within  2%  in  the  magnitude  and  2°  in  phase. 

To  verify  that  the  size  of  the  aperture 
does  not  significantly  affect  the  results  or 
require  a  larger  ABC  spacing,  consider  the  data 
shown  in  Figure  7.  This  is  for  a  waveguide  of 
twice  the  dimensions  as  an  X-Band  waveguide 
(a=1.8"  and  b=0.8").  It  is  operated  at  10 
GHz.,  making  this  antenna  twice  as  large 
electrically  as  the  one  simulated  previously. 
For  this  case  the  first  6  modes  in  the  aperture 
are  compared.  The  ABC  was  again  rectangular 

and  was  spaced  7J2it  from  the  aperture  and  the 
edges  of  the  ground.  In  this  case  there  is  also 
excellent  agreement  between  the  two  sets  of 
data,  the  dominant  mode  amplitudes  are 
different  by  less  than  5%. 

As  a  final  test,  consider  the  data 
presented  in  Figure  8.  For  this  case,  there  are 
two  X  band  waveguide  fed  apertures  (same 
dimensions  as  was  used  for  Figures  5  and  6)  in 
the  ground  plane.  The  geometry  is  shown  in 
Figure  8a.  The  apertures  are  spaced  0.6" 
apart.  One  is  driven  and  the  other  is  connected 
to  a  matched  load.  Due  to  mutual  coupling 
betw'een  the  two  aperture  antennas,  there  is  a 
non-zero  field  in  the  aperture  of  the  undriven 
guide.  The  modal  amplitudes  of  the  fields  in 
both  apertures  are  shown  in  Figure  8.  The 

ground  plane  again  extends  X/4  beyond  the 

apertures  and  the  ABC  spacing  is  7J2n. 

Conclusion 

From  the  data  presented  here  it  is  concluded 
that  the  new  second  order  ABCs  can  be  used 
very  effectively  with  the  finite  element  method 
to  model  aperture  radiation  problems.  The  new 
ABCs  can  be  shaped  to  conform  to  the  problem 

geometry  and  can  be  placed  as  close  as  X/2k 


from  the  aperture  for  accurate  simulation. 
Since  the  finite  element  method  is  well  suited  to 
modeling  real-life  structures  such  as  finite 
flanges  and  variations  in  enclosure  shapes,  it 
provides  a  new,  powerful  tool  to  analyze 
radiation  problems  efficiently. 

It  should  be  noted  that  these  results  are 
valid  for  the  types  of  geometries  presented 
here.  It  is  possible  that  a  complicated  radiating 
geometry  can  be  found  where  the  reactive  near 
zone  ends  at  a  distance  that  is  greater  than 

7J2n.  In  that  case,  the  ABC  should  be  spaced 
at  least  that  distance  from  the  geometry  for  best 
results. 
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Figure  1  A  tetrahedral  showing  the  electric  field  components  of  the  tangential  vector  finite  element. 


top  view  side  view 

Figure  2  FEM  simulation  model  of  the  unflanged  rectangular  waveguide  antenna. 
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2  (cm) 


Figure  3  Magnitude  of  the  near  zone  electric  field  along  the  center  axis  of  the  aperture.  The  values 
are  in  dBV/M. 

y 


infinite  ground  plane  (z=0) 

Figure  4  The  dimensions  and  coordinate  system  for  the  flanged  waveguide  antenna  used  in  the 
MOM  simulations.  For  the  X  band  waveguide  used  here:  a=0.9"  and  b=0.4". 
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Figure  5  Magnitude  of  the  modal  amplitudes  for  a  flanged  X  band  waveguide  antenna. 
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Figure  6  Phase  of  the  modal  amplitudes  for  a  flanged  X  band  waveguide. 
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Figure  7  Magnitude  of  the  modal  amplitudes  for  a  flanged  waveguide  antenna  that  is  twice  as  large 
as  a  standard  X  band  waveguide. 
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Figure  8  (a)  Two  aperture  geometry,  (b)  Magnitude  of  the  modal  amplitudes  for  the  case  of  a  X 
band  waveguides  with  d=0.6".  The  data  on  the  left  is  for  the  driven  guide,  and  the  data  on  the  right 
is  for  the  undriven  guide. 
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ABSTRACT.  The  structure  and  the  properties  of  the  FE¬ 
MAX  package  are  discussed.  The  FEMAX  finite-element 
package  is  an  efficient  and  highly  accurate  package 
especially  designed  for  computing  three-dimensional 
transient  as  well  as  time-harmonic  electromagnetic  fields 
in  arbitrarily  inhomogeneous j  (an)isotropic  media.  The 
most  unique  features  of  the  FEMAX  package  are  that 
1)  the  electric  field  strength  is  computed  directly,  i.e. 
without  the  intermediate  use  of  (vector)  potentials,  2) 
when  inhomogeneities  are  encountered  in  the  domain 
of  computation,  the  package  automatically  chooses  edge 
elements  to  ensure  that  all  local  continuity  conditions  can 
be  met,  nodal  elements  are  used  elsewhere,  and  3)  that 
the  electromagnetic  compatibility  relations  are  taken  into 
account  in  the  formulation  of  the  finite- element  method 
used,  thus  avoiding  spurious  solutions.  These  features 
are  included  in  the  package  in  such  a  way  that  optimum 
results  are  obtained  both  in  regard  to  computational 
efficiency  (storage  and  time)  and  in  regard  to  the  desired 
accuracy. 

1  INTRODUCTION 

When  solving  the  electromagnetic  field  equations  in  a 
three-dimensional  domain  containing  inhomogeneous  me¬ 
dia  by  using  finite  elements  some  important  issues  should 
be  dealt  with  that  are  not  often  encountered  in  finite- 
element  methods  for  other  sets  of  partial  differential  equa¬ 
tions. 

The  first  important  difference  between  the  electromag¬ 
netic  field  equations  and  other  sets  of  partial  differential 
equations  is  that  the  field  quantities  in  the  former  (the 
electric  and  the  magnetic  field  strength)  show  jump  dis¬ 
continuities  at  interfaces  between  different  media.  This 
makes  it  necessary  to  use  a  computational  technique  that 
accounts  for  the  continuity  conditions  of  the  fields  across 
interfaces  where  the  constitutive  coefficients  jump. 

The  second  issue  to  be  resolved  lies  in  the  fact  that 
finite-elements  methods  for  electromagnetic  field  prob¬ 
lems  often  yield  unwanted  (spurious)  solutions.  Those 


spurious  solutions  are  due  to  the  fact  that  some  of  the 
properties  of  the  electromagnetic  field,  the  electromag¬ 
netic  compatibility  relations  [1,  2],  are  not  represented 
properly  in  the  finite-element  method  used. 

1.1  Interfaces 

A  popular  method  to  simplify  the  modeling  of  the 
continuity  conditions  at  interfaces  is  to  introduce  electric 
and/or  magnetic  (vector)  potentials  [3]  that  are  con¬ 
tinuous  along  those  interfaces.  When  using  potentials 
in  a  numerical  method,  the  finite-element  method  for 
instance,  they  have  the  disadvantage  of  yielding  rela¬ 
tively  inaccurate  results  for  physical  quantities  like  the 
electric  or  magnetic  field  strength.  This  is  caused  by 
the  fact  that  the  latter  quantities  can  only  be  obtained 
by  performing  differentiations  on  the  numerical  results 
for  the  potentials  which  causes  the  loss  of  one  order  of 
accuracy.  Obviously,  much  better  convergence  properties 
are  obtained  when  numerical  differentiations  are  avoided 
by  formulating  the  problem  directly  in  terms  of  the 
electric  and/or  the  magnetic  field  strengths.  When 
doing  so,  the  continuity  conditions  along  interfaces  can, 
in  principle,  be  dealt  with  by  subdividing  the  problem 
space  into  a  number  of  subdomains  over  which  the 
constitutive  coefficients  are  continuous  functions  of  the 
spatial  coordinates.  The  boundary  conditions  at  the 
interfaces  between  those  subdomains  are  imposed  sepa¬ 
rately  [4,  5].  This  technique,  however,  is  very  difficult 
to  implement  for  arbitrary  configurations  and  yields 
conflicting  conditions  at  nodes  where  the  vector  normal 
to  the  interfaces  is  not  unique.  When  the  electromagnetic 
field  is  computed  in  terms  of  the  electric  and/or  the 
magnetic  field  strength,  the  problem  of  modeling  fields 
in  inhomogeneous  media  can  be  solved  at  the  element 
level  by  using  edge  elements.  Edge  elements,  however, 
are  known  to  be  less  efficient  than  the  commonly  used 
nodal  elements  [6].  In  [7]  a  method  is  described  that 
virtually  eliminates  the  computational  disadvantage  of 
edge  elements  [6]  by  using  them  only  for  modeling  the 


field  along  discontinuities  and  near  reentrant  corners, 
and  by  using  nodal  elements  everywhere  else.  With 
this  ’’mixed”  method  the  corresponding  program  decides 
locally,  for  each  combination  of  two  adjacent  subdomains 
in  the  discretized  domain  of  computation,  which  type  of 
element  will  be  used.  For  regions  in  which  the  properties 
of  the  media  are  continuous  functions  of  the  spatial 
variables  and  for  regions  with  ’’weak”  discontinuities 
(i.e.  discontinuities  along  which  ignoring  the  jump  in  the 
normal  component  of  the  field  strength  would  not  yield 
unacceptably  large  errors)  in  electromagnetic  properties 
between  subdomains  it  will  use  nodal  elements,  for 
’’large”  differences  it  will  use  edge  elements.  Thus  the 
method  of  modeling  the  field  is  automatically  adapted  to 
the  problem  at  hand  and  all  ’’strong”  discontinuities  are 
automatically  accounted  for.  The  degree  of  discontinuity 
above  which  edge  elements  are  used  is  user- defined  and 
has  to  be  chosen  in  accordance  with  the  final  accuracy 
aimed  at.  In  [8]  a  finite  element  code  for  time-harmonic 
electromagnetic  fields  was  described  that  was  designed 
according  to  the  method  described  above. 

1.2  Spurious  solutions 

Spurious  solutions  can  be  eliminated  by  taking  into 
account  all  electromagnetic  compatibility  in  the  formu¬ 
lation  of  the  finite-element  method  [1].  In  this  way  the 
method  used  will  accurately  model  the  electromagnetic 
field  equations,  together  with  all  its  relevant  properties, 
and  spurious  solutions  are  avoided. 

In  the  present  paper  the  version  of  the  FEMAX  finite- 
element  code  is  described  that  was  especially  designed 
for  computing  transient  electric  fields  in  three  spatial 
dimensions.  It  was  developed  on  the  basis  of  the  above 
ideas. 

2  NODAL  ELEMENTS  AND  EDGE  ELEMENTS 

As  was  explained  above,  the  FEMAX  package  uses  a 
combination  of  nodal  and  edge  elements.  For  topolog¬ 
ical  reasons  [9]  the  geometrical  domain  V,  in  which  the 
finite-element  method  is  applied,  is  subdivided  into  tetra- 
hedra  that  together  span  the  polyhedron  approximating 
D.  Consequently,  the  nodal  and  edge  elements  to  be  used 
will  be  defined  on  tetrahedra.  The  position  vectors  of  the 
vertices  of  a  particular  tetrahedron  T  are  {ro,  ,  r2,  rs}, 
the  outwardly  directed  vectorial  areas  of  the  faces  of  T 
are  {Ao,Ai,A2,A3}  and  the  volume  of  T  is  denoted  by 
V.  Let  Tb  be  the  position  vector  of  the  barycenter  of 
T,  then  the  linear  functions  ^i(r)  that  equal  unity  when 
r  =  xi\  +  yi2  +  zi^  =  Ti{i  =  0, 1, 2, 3)  and  are  zero  at  the 


remaining  three  vertices  of  T  can  be  written  as 

(j)i{r)  =  1/4  -  (r  -  Tb)  •  AilW.  (1) 

For  obtaining  consistency  in  the  local  degree  of  approx¬ 
imation  we  use  consistently  linear  expansion  functions. 
Therefore,  the  local  nodal  expansion  functions 
in  T  are  taken  as 

(^)  =  =  0, 1,2,3  and  j  =  1,2,3,  (2) 

and  the  local  edge  expansion  functions  W\  j{r)  in  T  are 
taken  as 

(^)  =  -Mr)aijAj/3V,  i,j  =  0, 1,2,3,  i  ^  j,  (3) 

where  aij  =  \ri  —  rj  \  is  used  for  dimensioning  and  scaling 
purposes.  When  r  £  T,  the  electric  field  strength  E  is 
locally  expanded  in  terms  of  a  combination  of  nodal  and 
edge  expansion  functions  through 

(4) 

i=Q  j 

where  eij{t)  are  the  unknown,  time-dependent,  field  ex¬ 
pansion  coefficients  and  where  the  limits  of  j  depend  on 
the  type  of  expansion  used  (see  (2)  and  (3)).  Note  that  we 
have,  in  each  tetrahedron,  the  possibility  to  use  arbitrary 
combinations  of  nodal  and  edge  expansion  functions.  The 
actual  choice  between  them  is  made  by  the  program,  it 
depends  on  the  degree  of  inhomogeneity  in  the  material 
properties  in  the  immediate  vicinity  of  T.  Nodal  expan¬ 
sion  coefficients  are  the  Cartesian  components  of  the  field 
at  the  vertices  of  the  tetrahedra  (the  nodes  of  the  mesh). 
Edge  expansion  coefficients  are  the  oriented  projections  of 
the  field  at  the  edges  taken  near  the  ends  of  those  edges. 
We  note  that  (4)  always  contains  12  terms.  For  a  dis¬ 
cussion  of  the  properties  of  local  expansion  functions  the 
reader  is  referred  to  [6]. 

We  now  introduce  global  expansion  functions  as  the 
sum  of  all  local  expansion  functions  that  are  related  to 
the  same  expansion  coefficient.  The  support  of  the  global 
nodal  expansion  functions  is  the  simplicial  star  [9]  of  the 
node  to  which  they  are  related  (the  simplicial  star  of  a 
node  is  the  union  of  all  tetrahedra  that  have  that  node 
in  common).  The  support  of  the  global  edge  expansion 
functions  is  the  simplicial  star  of  the  edge  to  which 
they  are  related  (the  simplicial  star  of  an  edge  is  the 
union  of  all  tetrahedra  that  have  that  edge  in  common). 
Since  the  support  of  global  nodal  expansion  functions  is 
much  larger  than  the  support  of  global  edge  expansion 
functions,  the  former  are  computationally  more  efficient 
than  the  latter  and  the  former  should  be  used  whenever 
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the  use  of  the  latter  is  not  necessary  because  of  the  local 
degree  of  inhomogeneity.  For  an  illustration  the  reader  is 
referred  to  Fig,  2  in  [8].  By  using  the  proper  combination 
of  nodal  and  edge  expansion  functions,  optimum  results 
are  obtained  both  in  regard  to  accuracy  and  in  regard  to 
computational  efficiency  [6,  8]. 

3  THE  FINITE-ELEMENT  FORMULATION 


3.1  The  wave  equation 

Eliminating  the  magnetic  field  strength  H  from  Max¬ 
well’s  equations  we  obtain 


die  -  E  +  dt(T  ■  E  +  V  X  -  V  X  E)  = 

-  V  X  (/x-i  •  (5) 

where  =  J™P(r,<)  and  denote 

known  impressed  sources  of  electrical  and  magnetic  cur¬ 
rent  that  are  defined  inside  the  domain  of  computation 
only,  and  where  e  —  €(r),  a  =  cr(r)  and  /z  =  |x(r) 
denote  the  permittivity,  the  conductivity  and  the  perme¬ 
ability  tensor,  respectively.  The  wave  equation  is  sup¬ 
plemented  with  initial  conditions  E(r,io)  and  i/(r,^o)- 
In  the  package  the  domain  of  computation  is  subdivided 
in  tetrahedra.  Now,  substituting  the  expansion  (4)  for 
the  electric  field  strength  in  (5),  a  system  of  equations 
in  the  expansion  coefficients  is  obtained  by  applying  the 
method  of  weighted  residuals.  The  set  of  weighting  func¬ 
tions  that  is  used  is  the  same  as  the  set  of 

expansion  functions.  Using  an  integration  by  parts  and 
adding  the  resulting  equations  over  all  tetrahedra,  we  ob¬ 
tain  a  system  of  coupled  ordinary  differential  equations 
for  {cij{t)}  that  can  be  written  as 

[  Wp,,-€-WijdV 

Jv 


Figure  1:  The  domain  of  computation  V. 


■  Vp,q, 


(6) 


where  dVn,  see  Fig.  1,  denotes  the  part  of  outer  bound¬ 
ary  of  the  domain  of  computation  V  on  which  the  tan¬ 
gential  components  of  the  external  magnetic  field  strength 
t)  are  prescribed  and  where  n  denotes  the 
unit  vector  along  the  outward  normal  to  dV,  The  tan¬ 
gential  components  of  the  external  electric  field  stength 
—  E^^^{r,t)  are  prescribed  on  where  dV  = 

dVE  U  dVu  with  OVe  H  dDu  —  0.  External  fields  are 
defined  outside  the  domain  of  computation  only.  For  de¬ 
riving  (6)  we  have  used  the  continuity  of  the  tangential 
components  of  the  magnetic  field  strength  over  all  inter¬ 
nal  interfaces. 


3.2  The  compatibility  relations 

For  eliminating  all  spurious  solutions  (6)  must  be  aug¬ 
mented  with  the  relevant  electromagnetic  compatibility 
relations  [1].  Compatibility  relations  are  properties  of  a 
field  that  are  direct  consequences  of  the  field  equations 
and  that  should  be  satisfied  to  allow  the  equations  to 
have  a  solution.  They  were  first  introduced  by  Love  [10] 
for  elastodynamics  and  are  also  known  in  fluid  dynamics 
[11].  For  the  electric  field  strength,  which  is  the  fun¬ 
damental  unknown  in  FEMAX,  these  compatibility  rela¬ 
tions  are  the  volume  divergence  condition 


+  /  Wp,,  ■  (T  •  WijdV 

ij 

+  /  (V  X  Wp,,)  ■  /i-i  •  (V  X  Wij)dV  = 

id 

[  Wp^,  ■  in  X  H^^^)dA 

Jv 


V  •  -  E  -f-  cr  ^  E)  =  -  V  •  (7) 

the  surface  divergence  condition  at  interfaces 

u  •  {die  '  E  a  ’  E)  1/  ‘  continuous  across  X,  (8) 

where  u  denotes  the  unit  vector  along  the  normal  to  the 
interface  X,  and  the  condition  applying  to  the  normal 
component  of  the  electric  flux  density  at  the  outer  bound¬ 
ary  dVn 
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n-{die^E^a^E)  =  7i^{VxW^^-J^^^)  ondVu^  (9) 

Each  of  these  equations  is  added  to  (6)  in  a  weighted 
sense.  In  this  way  the  wave  equation,  together  with  its 
properties,  is  modelled  in  a  weighted  sense  and  spurious 
solutions  are  avoided. 

Adding  the  weighted  form  of  (7)- (9)  to  (6),  a  system 
of  coupled  ordinary  differential  equations  for  {eij{t)] 
is  obtained  that  is  solved  in  a  standard  manner  using 
either  single-  [12]  or  two-step  [13]  integration  methods. 
In  practice  the  use  of  two-step  methods  has  turned  out 
to  be  preferable. 

4  TIME-HARMONIC  FIELDS  IN  THE  TIME 
DOMAIN 

The  FEMAX  code  allows  time-harmonic  fields  to  be  com¬ 
puted  in  the  time  domain.  For  maximum  efficiency  a  tran¬ 
sient  is  used  with  a  transient  period  0  <  i  <  ttr-  During 
this  period  the  known  right-hand  side  vector,  representing 
the  time  harmonic  excitation,  is  multiplied  by  a  continu¬ 
ous  function  /tr(i)  that  generates  a  smooth  transient  from 
zero  to  steady  state.  For  fu  we  use  a  function  of  the  type 

=  0,  -00  <t  <0, 

ftr  =  (2-sin(j^f))sm(j^f),  0  <  f  <  ftr,  (10) 

—  1,  ttr  t  00, 

other  transients  being  available.  Using  /tr  we  obtain, 
at  t  =  ttrj  an  approximation  of  the  time-harmonic  so¬ 
lution  with  a  relatively  small  error  term  and  steady  state 
is  achieved  very  efficiently.  In  most  cases  a  value  of  ^tr  in 
the  range  5T  <  ttr  <  lOT,  where  T  denotes  the  period  in 
time  of  the  time*harmonic  sources,  turns  out  to  yield  an 
optimum  computational  efficiency.  Obviously  the  opti¬ 
mum  choice  for  ttr  depends  both  on  the  problem  at  hand 
(larger  domains  of  computation  requiring  longer  transient 
times)  and  on  the  accuracy  requirements.  In  general  it 
can  be  said  that  ttr  should  preferably  be  chosen  such  that 
the  solution  at  t  =  ttr  is  already  accurate  enough  to  be 
used  as  a  steady  state  solution  and  that  no  subsequent 
time  stepping  is  required  to  obtain  steady  state.  For  fur¬ 
ther  comments  on  this  procedure  the  reader  is  referred  to 
[14]. 

Computing  time-harmonic  fields  in  the  time  domain 
is  computationally  very  attractive  as  compared  with 
using  frequency  domain  methods  because  of  the  fact  it 
is  usually  very  difficult  to  solve  the  system  of  equations 
obtained  using  the  latter  method  which  tends  to  generate 


ill  conditioned  system  matrices  that  usually  are  not 
positive  definite.  Experimentally  we  have  found  that  for 
’’large”  problems  (>10.000  unknowns)  the  time  domain 
approach  to  solving  time-harmonic  problems  is  always 
more  efficient  and  more  accurate  than  the  time-harmonic 
approach.  The  advantage  of  the  time-domain,  transient, 
approach  over  the  time-harmonic  approach  increases 
with  the  number  of  unknowns  and  with  the  temporal  fre¬ 
quency  (for  a  constant  size  of  the  domain  of  computation). 

5  ABSORBING  BOUNDARY  CONDITIONS 

The  FEMAX  package  contains  simple  local  absorbing 
boundary  conditions  that  model  the  infinite  homoge¬ 
neous  surroundings  of  the  domain  of  computation.  These 
boundary  conditions  assume  the  field  locally  to  consist 
of  a  plane  wave  travelling  in  a  specific  direction,  usually 
in  the  direction  that  is  normal  to  the  outer  boundary, 
out  of  the  domain  of  computation.  Absorbing  boundary 
conditions  of  this  type  can  be  used  for  problems  in  which 
all  sources  of  the  electromagnetic  field  are  located  inside 
the  domain  of  computation. 

5.1  Inhomogeneous  absorbing  boundary  condi¬ 
tions 

In  many  practical  problems,  the  field  is  excited  by  sources 
that  are  located  outside  the  domain  of  computation.  As 
an  example  one  might  think  of  the  scattering  of  a  plane 
wave  by  an  obstacle.  In  cases  like  this  we  have  an  inci¬ 
dent  field  which  is  the  field  in  the  absence 

of  the  scatterer,  a  total  field  {E,  JT},  and  a  scattered 
field  which  is  defined  as  the  difference  be¬ 

tween  the  total  field  and  the  incident  field.  Absorbing 
boundary  conditions  apply  to  scattered  fields,  they  can 
be  augmented  to  take  into  account  the  incident  field  in 
the  following  manner.  Since  absorbing  boundary  condi¬ 
tions  are  linear  expressions  in  terms  of  the  scattered  field 
they  can  formally  be  written  as  a  linear  operator  of  the 
type 


|j(^scat  j^scat)  ^  0. 

(11) 

Adding  the  trivial  identity 

(12) 

to  (11)  we  obtain 

i?(jE;,  jf)  =  «■■"*=), 

(13) 

which  relation  is  an  inhomogeneous  absorbing  boundary 
condition  that  can  be  applied  to  the  total  field  inside  the 
domain  of  computation  and  that  takes  into  account  the 
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incident  field  as  a  contribution  to  the  known  excitation 
vector.  Note  that  (13)  applies  to  all  absorbing  boundary 
conditions  that  can  be  written  as  (11),  from  very  simple 
local  conditions  to  highly  sophisticated  global  ones. 

6  THE  PROPERTIES  OF  THE  FEMAX  PACK¬ 
AGE 

The  FEMAX  package  for  computing  time-domain  elec¬ 
tromagnetic  fields  in  three  spatial  dimensions  has  the  fol¬ 
lowing  properties  and/or  options: 

•  The  electric  field  strength  is  used  as  the  fundamental 
unknown  quantity  that  will  be  approximated  using 
consistently  linear  expansion  functions. 

•  Because  of  including  all  relevant  electromagnetic 
compatibility  relations  in  the  formulation,  the  results 
are  free  of  spurious  solutions. 

•  In  the  interior  of  the  domain  of  computation,  the  user 
can  specify  arbitrary  distributions  of  electromagnetic 
medium  properties  (permittivity,  conductivity  and 
permeability),  the  media  may  be  anisotropic.  It  is 
also  possible  to  specify  specific  density  distributions 
for  use  in  SAR  computations. 

•  The  package  uses  edge  elements  along  discontinu¬ 
ities  and  nodal  elements  in  regions  with  continuously 
varying  medium  properties.  The  contrast  in  prop¬ 
erties  of  the  media  in  two  adjacent  domains  below 
which  the  medium  is  treated  as  having  continuously 
varying  properties  is  specified  by  the  user. 

•  For  a  given  distribution  of  medium  properties  and 
a  given  finite-element  mesh,  which  has  to  be  chosen 
fine  enough  to  satisfy  the  accuracy  requirements,  the 
package  automatically  determines  the  numerically 
optimum  distribution  of  the  expansion  functions  that 
allow  the  solution  to  satisfy  the  user-defined  local  ac¬ 
curacy  requirements. 

•  In  the  interior  of  the  domain  of  computation,  the 
user  can  prescribe  arbitrary  initial  value  distribu¬ 
tions  of  the  field  strength  as  well  as  arbitrary  time- 
dependent  volume  source  distributions  of  imposed 
electric  and/or  magnetic  current. 

•  The  distributions  of  the  medium  properties  as  well 
as  the  volume  source  distributions  can  be  specified 
either  by  supplying  user-written  subroutines  for  gen¬ 
erating  those  quantities  (in  FORTRAN- 77)  or,  much 
simpler,  by  using  members  of  the  collection  of  stan¬ 
dard  subdomains  with  standard  distributions  that 
are  available  in  the  package. 


•  At  the  outer  boundary  of  the  domain  of  computa¬ 
tion,  the  user  can  prescribe  arbitrary  distributions 
of  external  electromagnetic  fields.  Local  impedance 
boundary  conditions  for  modeling  lossy  boundaries 
and  (in)homogeneous  absorbing  boundary  conditions 
for  modeling  radiation  into  an  unbounded  homoge¬ 
neous  medium  surrounding  the  domain  of  computa¬ 
tion  are  also  available. 

•  Efiicient  sparse-matrix  methods  and  different  types 
of  preconditioning  are  available  for  solving  the  sys¬ 
tem  of  algebraic  equations. 

•  Field  values  and  derived  quantities  can  be  computed 
at  arbitrary  collections  of  user-specified  points  in  the 
configuration  and  as  a  function  of  time. 

•  FEMAX  uses  the  SEPRAN  finite-element  package 
[18]  for  a  number  of  general  finite-element  tasks  like 
the  generation  of  the  mesh  and  the  assembly  of  the 
system  matrix  from  the  FEMAX  elements. 

A  number  of  programs  is  available  for  post-processing 
data  prepared  in  FEMAXT.  They  can  be  used  inter¬ 
actively  for  previewing  and  for  making  contour  plots 
(FEMAXPT  and  FEMAXPF)  of  field  values,  and  of 
quantities  that  can  be  expressed  in  terms  of  the  field 
values  and  material  properties.  FEMAXH  is  used  for 
making  plots  of  the  evolution  in  time  of  fields  values 
at  specified  locations.  Use  is  made  of  either  the  NAG 
Graphics  Library  [15]  or  MATLAB  [16]. 

7  NUMERICAL  RESULTS  AND  COMPUTA¬ 
TIONAL  REQUIREMENTS 

Early  results  obtained  by  using  the  FEMAX  package  have 
been  presented  in  [14,  17].  The  results  presented  in  those 
papers  mainly  served  to  illustrate  the  accuracy  and  the 
convergence  properties  of  the  package  on  a  simple  model 
problem.  As  a  more  practical  example  we  shall  now  com¬ 
pute  the  specific  absorption  rate  of  a  field  in  a  mathemat¬ 
ical  model  of  a  human  head  (a  numerical  phantom).  The 
field  is  generated  by  a  dipole  antenna  that  is  located  near 
this  phantom.  A  quarter  part  of  the  configuration  is  de¬ 
picted  in  Fig.  2.  The  frequency  is  taken  as  /=900MHz, 
which  is  inside  the  frequency  range  used  for  the  Euro¬ 
pean  GSM  (Global  System  for  Mobile  Communication) 
network  range.  The  phantom  has  a  cubic  shape,  it  con¬ 
sists  of  a  cubic  region  -0.095m.  <  x^y^z  <  0.095m  filled 
with  brain  tissue  (£r=43,  a  =0.83S/m,  p  =1050kg/m^). 
The  region  outside  the  above  region,  but  inside  the  cube 
-0.1m  <  x,y,z  <  0.1m  has  properties  for  modelling  the 
skull  (5r=17,  (T  =0.25S/m,  p  =1200kg/m3).  Outside  this 
region  the  properties  of  vacuum  are  assumed. 
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Figure  2:  The  numerical  phantom 

The  antenna  is  modelled  by  a  A/2  dipole  of  length 
0.166m,  having  its  center  at  {x,y,z)  =  (0.115,0,0)  and 
its  orientation  parallel  to  the  z-axis. 

The  problem  was  solved  under  the  condition  that  edge 
elements  are  used  when  the  relative  contrast  Cr  in  the  nu¬ 
merical  values  of  cr(r)-fjcc;er(T’)eo  in  two  adjacent  tetrahe- 
dra  exceeds  0.1  (with  this  choice,  edge  elements  are  used 
along  all  interfaces  between  the  different  domains),  nodal 
elements  are  used  otherwise.  By  using  the  symmetry  of 
the  configuration  with  respect  to  the  planes  y  =  0  and 

=  0  we  can  solve  the  problem  for  positive  y  and  z  values 
only,  using  the  condition  E^[r^t)  E{r^t)  —  0  (the 

subscript  T  refers  to  the  tangential  components)  at  the 
plane  z  =  0  and  the  condition  Ht  (r,  f)  =  ^2  x  i/(r ,  i)  =  0 
at  the  plane  y  =  0.  The  domain  of  computation  —0.11  < 
X  <  0.2,0  <  y  <  0.16,0  <  z  <  0.16  is  surrounded  by 
a  mathematical  boundary  on  which  absorbing  boundary 
conditions  are  applied.  For  the  problem  at  hand  the  solu¬ 
tion  is  expected  to  have  a  ’’singularity”  in  the  immediate 
vicinity  of  the  dipole  antenna.  Therefore  we  have  used  a 
non-uniform  mesh  that  has  its  highest  density  near  this 
antenna.  A  side  view  of  this  mesh,  consisting  of  144000 
tetrahedra,  is  given  in  Fig.  3.  Note  that  the  mesh  is 
such  that  the  interfaces  between  different  media  coincide 
with  interfaces  between  tetrahedra,  in  this  way  we  have 
avoided  staircasing  in  the  modeling  of  the  material  prop¬ 
erties. 

For  this  configuration  the  Specific  Absorption  Rate 


(SAR  in  W/kg)  is  required.  In  Fig.  4  the  normalized 
SAR  distribution  is  depicted  for  a  cross-section  in  the 
plane  y=0.  The  normalization  was  carried  out  such  that 
the  maximum  SAR  value  anywhere  in  the  configuration 
was  set  to  OdB.  In  Fig.  5  the  SAR  distribution  is  de¬ 
picted  along  the  line  (-0.1  <  x  <  0.1, y  =  0,z  =  0),  the 
same  normalization  was  used.  The  almost  linear  decay  of 
the  SAR  (dB)  values  with  decreasing  x-values  reflects  the 
exponential  decay  of  the  field  in  the  highly  lossy  tissues. 
The  ’’roughness”  of  the  result  in  Fig.  5  is  due  to  the  fact 
that  a  relatively  coarse  mesh  had  to  be  used  away  from 
the  antenna.  For  x  <  0.05  about  5  elements  per  wave¬ 
length  were  used,  near  the  antenna  a  much  denser  mesh 
was  used  and,  consequently,  much  smoother  results  were 
obtained  for  a;  >0.05. 


Fig.  3.  The  mesh  used 
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The  computational  effort  for  solving  a  problem  follows 
from  the  effort  required  for  generating  the  system  matri¬ 
ces,  the  effort  for  solving  the  equations  each  time  step 
and  the  number  of  time  steps.  Since  the  matrices  are 
stored  in  compact  mode  (ignoring  all  zero’s)  the  solution 
of  the  system  of  equations  can  be  obtained  only  by  using 
iterative  methods  because  of  which  the  solution  time  de¬ 
pends  on  the  number  of  iterations  required.  The  actual 
number  of  iterations  for  solving  a  problem  iteratively  each 
time  step  is  relatively  small  and  depends  on  the  type  of 
preconditioning  available  for  the  problem  at  hand.  Typ¬ 
ically  2-20  iterations  are  required  for  solving  a  system 
of  equations  with  100,000  unknowns.  The  computational 
effort  increases  only  slightly  faster  than  linear  with  the 
number  of  unknowns  and  linear  with  the  number  of  time 
steps  required.  For  large  problems  with  many  time  steps 
the  computational  effort  for  generating  the  system  matri¬ 
ces  is  negligible  in  comparison  to  the  computational  effort 
required  for  the  integration  of  the  system  of  differential 
equations  along  the  time  axis. 


Fig.  4.  The  normalized  SAR  distribution  in  the  plane 
7/  =  0  in  dB. 


For  the  problem  at  hand  the  total  number  of  un¬ 
knowns  amounted  to  114085,  3082  of  them  being 
prescribed  through  essential  boundary  conditions.  Each 
of  the  three  matrices  contained  4116764  non-zero  ele¬ 
ments.  The  time-harmonic  problem  was  solved  with  a 
transient  of  10  periods  in  time,  with  20  time  steps  on 
each  period  a  total  of  1055  iterations  were  required  to 
solve  the  problem.  All  computations  were  carried  out  on 
a  HP  9000-735/125  workstation,  with  400Mbyte  main 
memory  space.  The  total  computation  time  for  obtaining 
the  solution  amounted  to  about  100  min. 


8  THE  PACKAGE 

The  FEMAX  finite-element  package  consists  of  approxi¬ 
mately  250  subroutines.  It  uses  the  SEPRAN  [18]  finite 
element  package  as  a  general  finite-element  background 
package  that  carries  out  tasks  like  the  generation  of  the 
mesh,  the  assembly  of  the  matrices  and  vectors,  the  so¬ 
lution  of  the  systems  of  algebraic  equations  and  the  data 
management.  Either  the  NAG  Graphics  Library  [15] 
or  MATLAB  [16]  is  used  for  postprocessing.  The  stor¬ 
age  requirements  for  the  executable  are  6.5Mbytes.  FE¬ 
MAX  runs  on  both  VAX- VMS  and  UNIX  platforms.  The 
FEMAX  source  code,  together  with  an  extensive  User’s 
Guide,  are  commercially  available.  The  SEPRAN  source 
code  should  be  acquired  separately. 

The  time-harmonic  version  of  the  FEMAX  packages 
[8]  is  also  commercially  available.  A  FEMAX  package  for 
static  and  stationary  electric  or  magnetic  fields  based  on 
the  theory  presented  in  [2]  is  presently  under  development 
[19]. 
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Fig.  5.  The  normalized  SAR  distribution  along  the 
line  t/  =  0,z  =  0  in  dB. 

9  CONCLUSION 

The  FEMAX  finite-element  package  was  described.  This 
package  was  especially  designed  for  computing  transient 
as  well  as  time-harmonic  electromagnetic  fields  in  three- 
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dimensional  configurations  containing  inhomogeneous 
and  (an)isotropic  media.  Because  of  formulating  the 
problem  in  terms  of  the  electric  field  strength  directly,  all 
problems  that  axe  inherent  to  potential  formulations  are 
avoided.  The  modeling  of  the  conditions  along  internal 
interfaces  in  the  configuration  is  performed  automatically 
by  the  package,  and  the  user  can  infiuence  this  process 
by  prescribing  the  minimum  degree  of  discontinuity  Cr 
for  which  the  discontinuity  in  the  medium  properties 
should  be  modelled  to  obtain  the  required  accuracy. 
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ABSTRACT  A  Moment  Method  (MM)  technique  for  the 
analysis  of  wire  antennas  on  board  resonant-sized  bodies 
modelled  with  parametric  surfaces  is  presented.  The 
approach  may  have  useful  applications  for  the  study  of  the 
behaviour  of  wire  antennas  on  board  complex  conducting 
structures  like  aircrafts.  The  current  is  represented  by 
curved  rooftop  functions  on  the  body,  piecewise  linear 
functions  on  the  wires  and  a  new  junction  function  in  the 
attachment  region  between  the  body  and  the  antenna.  The 
bodies  are  precisely  defined  by  means  of  a  small  number  of 
NURBS  surfaces  (Non  Uniform  Rational  B-Splines)  and 
Bezier  patches  (BP).  In  addition  the  new  junction  function 
that  can  be  defined  over  any  BP  allows  the  antenna  to  be 
attached  to  any  part  of  the  body.  Radiation  patterns  and 
input  impedance  calculations  for  several  geometries  are 
presented  to  show  the  accuracy  of  the  method.  The  results 
are  successfully  validated  when  comparisons  with 
measurements  or  results  from  other  methods  are  carried  out. 

1  INTRODUCTION 

The  MM  [1]  is  one  of  the  most  popular  techniques  for  the 
analysis  of  electromagnetic  scattering  of  arbitrary  bodies  of 
small  or  resonant  size.  The  method  gives  good  results  when 
analyzing  wire  antennas  and  can  be  used  to  obtain  the 
scattering  field  by  surfaces  represented  by  wire  grid  models 
[2].  However,  the  wire  grid  modelling  does  not  work 
properly  when  near-field  parameters  such  as  input 
impedances  or  current  distributions  are  calculated. 
Subsequently,  a  surface  representation  of  the  conducting 
bodies  must  be  introduced. 

A  representation  of  the  bodies  by  flat  triangular  or 
rectangular  patches  provides  excellent  results  of  near  and 
far-field  parameters  for  many  problems,  however  it  presents 
the  disadvantage  that  high  amounts  of  memory  and  CPU- 
times  are  needed  when  the  geometry  is  electrically  large. 
These  problems  can  be  substantially  reduced  by  defining  the 
MM  basis  functions  directly  over  the  models  of  the  bodies 
based  on  meshes  of  NURBS  surfaces  and  BP  [3]  because 
fewer  number  of  patches  are  needed  to  represent  the  body 
accurately. 

NURBS  surfaces  are  B-Spline  elements  formed  by  a  set  of 
patches  defined  by  polynomials,  [3].  The  coefficients  of 


these  polynomials  depend  on  merely  a  few  control  points 
(see  Fig.l).  A  body  that  is  quite  complex,  for  instance  a 
complete  aircraft,  can  be  modelled  with  nearly  all  its  details 
by  only  a  few  hundred  NURBS.  Any  NURBS  can  be 
expressed  in  terms  of  BP.  NURBS  surfaces  are  more 
efficient  for  representing  and  storing  the  information 
necessary  for  describing  a  geometry  however  BP  are  more 
suitable  for  the  numerical  computation  of  parameters 
associated  with  the  local  behaviour  of  a  surface  like 
curvatures,  derivatives,  etc.  It  is  a  fast  and  easy  process  to 
obtain  the  BP  fi-om  a  NURBS  representation  by  applying  the 
Cox-de  Boor  transformation  algorithm  [4]. 


control  points 


Figure  1.  NURBS  surface  and  its  associated  control  points 


In  a  previous  paper  [5]  the  authors  presented  a  MM 
technique  using  NURBS  surfaces  and  BP  to  represent  the 
geometry  of  arbitrarily  shaped  scatterers.  Modelling  using 
that  method  implies  that  fewer  basis  functions  are  needed  to 
obtain  a  precise  representation  of  complex  bodies.  In  that 
approach  a  new  basis  function  associated  with  each 
boundary  line  between  pairs  of  adjacent  BP  was  introduced. 
These  basis  functions  can  be  considered  a  generalization  of 
the  planar  rooftop  functions  introduced  by  Glisson  [6]. 
Curved  ’’razor-blade”  functions  were  considered  as  testing 
functions.  This  method  was  successfully  validated  when 
RCS  values  of  several  objects  were  obtained  and  compared 
with  available  data  from  other  methods. 

The  purpose  of  this  paper  is  to  present  an  extension  of  the 
method  mentioned  above  in  order  to  analyze  antennas 
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mounted  over  conducting  surfaces.  By  establishing 
continuity  of  current  at  the  wire-surface  junction,  a  new 
special  basis  function  (attachment  or  junction  subdomain) 
has  been  introduced.  This  junction  subdomain  extends  over 
a  BP  of  the  model  and  a  segment  of  the  attached  wire 
antenna.  This  new  subdomain  presents  an  added  advantage 
compared  with  previous  works  that  can  be  attached  to  any 
of  the  BP  of  the  body  model.  As  the  BP  can  be  of  complex 
shape  it  is  possible  to  connect  the  antenna  to  any  part  of  the 
body.  Moreover  the  wire  attachment  itself  need  not  be  taken 
into  account  when  the  body  is  modelled.  Therefore  our 
formulation  derived  from  the  use  of  NURBS  and  BP  allows 
us  to  consider  more  complex  bodies  and  attachment  points 
than  in  previously  published  works,  [7-9]. 

This  paper  is  organized  as  follows:  Part  2  presents  the 
description  of  the  proposed  MM  scheme,  giving  an  overview 
of  the  basis  and  testing  functions  used,  the  junction 
subdomain,  the  electrical  fed  model  and  the  computation  of 
the  MM  matrix  terms;  Part  3  presents  the  radiation  patterns 
and  input  impedances  for  several  canonical  objects  including 
comparisons  with  results  obtained  by  other  methods  or 
measurements;  the  conclusions  are  outlined  in  Part  IV. 

2  FORMULATION 

The  electric  current  is  represented  by  three  kinds  of  basis 
functions  depending  on  the  part  of  the  geometry  considered: 
generalized  curved  rooftops  over  the  BP  that  describe  the 
body  surface,  piecewise  triangular  functions  on  the  wire 
antennas  and  junction  subdomains  in  the  wire-surface 
attachment  area. 

2.T  Generalized  curved  rooftops 

The  body  is  modelled  by  means  of  meshes  of  NURBS 
surfaces  that  are  subdivided  into  BP.  A  basis  function  is 
assigned  to  each  one  of  the  boundary  lines  between  pairs  of 
adjacent  BP;  each  basis  function  extends  only  over  the  pair 
of  patches  that  share  a  common  boundary  line.  A  constant 
charge  density  is  imposed  on  the  BP,  [5].  Curved  razor- 
blade  functions,  defined  over  the  isoparametric  lines  which 
join  the  centres  of  the  pair  of  BP  associated  with  each  basis 
function  have  been  considered  as  testing  functions. 

2.2  Wire  antennas 

A  thin  wire  approximation  has  been  used,  assuming  that  the 
radius  of  the  antennas  is  much  smaller  than  the  length. 
Although  the  wires  can  be  arbitrarily  curved,  a  piecewise 
linear  approximation  is  made  by  connecting  straight 
cylindrical  sections  as  can  be  seen  in  Fig.2.  Triangle  and 
pulse  functions  have  been  chosen  as  basis  and  testing 
functions,  respectively.  A  procedure  which  bears  a 


considerable  resemblance  to  that  proposed  in  [6]  has  been 
considered  for  the  numerical  treatment  of  the  MM  terms. 


Figure  2.  Curved  antennas  subdivided  into  a  set  of 
cylindrical  straight  segments 

2.3  Attachment  subdomain 

This  new  subdomain  extends  over  a  wire  segment  and  over 
a  BP,  as  can  be  seen  in  Fig.3.  The  wire  segment  is  assumed 
electrically  short  (e.g.  X/id  long)  and  normal  to  the  surface 
(the  other  wire  segments  can  be  arbitrarily  oriented).  The 
current  associated  with  the  junction  subdomain  is  a 
piecewise  linear  function  (semi-triangle)  on  the  wire  and  it 
is  assumed  that  the  current  flows  approximately  radially 
from  the  attachment  point  on  the  BP.  A  total  charge  of 
- 1  /y  0)  on  the  wire  (where  co  is  the  radian  frequency)  and 
a  charge  of  1  /y  co  on  the  patch  are  imposed.  The  charge 
continuity  equation  (/  =  -y  co  g)  tells  us  that  the  total  current 
of  the  junction  subdomain  is  1  Ampere.  The  charge  density 
is  assumed  to  be  constant  in  both  wire  segment  and  surface 
patch.  The  exact  shape  of  the  current  density  on  the  BP  of 
the  junction  subdomain  is  neither  defined  nor  computed. 
The  only  imposition  is  that  the  charge  density  should  be  a 
constant  function  with  a  total  charge  of  1  /y  cu .  As  we  will 
see  below  it  is  not  necessary  to  know  the  exact  shape  of  the 
current  on  the  BP  of  the  junction  subdomain. 


Figure  3.  The  junction  subdomain  formed  by  a  wire 
segment  and  a  Bezier  patch 
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Fig,4  shows  the  shape  of  the  testing  function  for  the 
junction  subdomain.  A  pulse  function  is  defined  as  a  testing 
function  solely  on  the  wire  part  of  the  junction  subdomain. 
No  testing  is  made  in  the  patch  surface  of  the  junction 
subdomain  because  this  patch  is  a  part  of  the  body 
geometry.  It  has  been  tested  by  the  weighting  functions 
associated  with  the  generalized  rooftop  functions  employed 
for  the  body,  [5]. 


Figure  4.  Testing  functions  for  the  junction  subdomain 


where  B  is  the  number  of  subdomains  over  the  patches  of 
the  geometry  modelled  with  NURBS  surfaces,  W  is  the 
number  of  subdomains  defined  on  all  the  wire  antennas  of 
the  problem  and  U  is  the  number  of  junction  subdomains. 
Then,  the  EFIE  for  this  problem  takes  the  form: 


«x£'V)=«xEse[-/?(r)] 


1  =  1 
w 


(2) 


i  =  1  /  =  1 


where  )  represents  the  impressed  electric  field,  n  is 

the  normal  vector  and  the  integro-differential  operator 

)]  is  defined  by: 


It  is  important  to  note  that  since  one  part  of  the  junction 
subdomain  is  one  patch  of  the  NURBS  meshing  (not  an 
additional  disc  or  a  flat  patch),  the  connection  of  an  antenna 
to  the  body  does  not  imply  any  restriction  or  modification 
in  the  NURBS  meshing.  The  sole  restriction  for  the 
attachment  points  is  that  they  should  be  positioned  in  or 
near  the  centre  of  a  BP.  This  is  not  an  important  limitation 
because  the  junction  patches  are  or  can  be  made  electrically 
small  (typically  a  body  is  modelled  using  four  to  six  patches 
per  wavelength  in  its  curved  areas  or  eight  to  ten  in  its  flat 
parts).  In  fact,  this  restriction  is  inherent  to  the  discretization 
process  of  the  problem  that  appears  when  using  a  numerical 
method  like  MM.  In  other  words  we  must  accept  a 
quantization  effect  for  the  possible  position  where  the  wire 
can  be  placed.  When  attaching  a  wire  antenna  to  a  patch 
model  of  a  body,  two  points  on  its  surface  are  considered 
different  only  if  they  do  not  belong  to  the  same  BP. 

2.4  MM  coupling  matrix 

The  objective  is  to  solve  the  electric  field  integral  equation 

(EFIE).  To  do  that,  the  current  density  J  is  expressed  as  a 
sum  of  the  current  density  on  each  subdomain  as: 


Equation  (2)  is  converted  into  an  equivalent  set  of  linear 
equations.  The  coefficients  of  this  system  form  the  MM 
impedance  matrix  that  can  be  grouped  in  nine  boxes: 


/z"“ 

z“* 

1 

z'^ 

^wb 

\z^  \ 

zbw 

2** 

where  u  denotes  a  junction  subdomain,  w  a  subdomain  on 
a  wire  and  b  a  subdomain  on  the  body. 

Since  the  objective  of  this  paper  is  to  introduce  the 
connection  between  wires  and  surfaces,  only  the  terms 
related  with  junction  subdomains  will  be  addressed  in  detail. 
The  matrix  terms  related  to  the  wire  antennas  have  been 
treated  by  other  authors,  for  instance  [7-9]  and  those  related 
to  the  subdomains  over  BP  were  introduced  in  [5]. 

If  we  consider  the  coupling  between  the  i  and  the  y-junction 

subdomains,  we  are  referring  to  one  of  the  terms  of  theZ““ 
box  of  the  MM  matrix.  This  term  is  obtained  by  evaluating 
the  expression: 


Ar)= E  ft  in  -  E  J7in  -  E  ftin 

i=l  i=l  i=l 


(1) 


z‘-{ij)  =  {Tt,^[j;]) 


(4) 
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where  7“  is  the  testing  function  associated  with  the  /- 
junction  subdomain  (see  Fig.4).  This  last  expression  can  be 
computed  as  a  sum  of  one  inductive  and  one  capacitive  term 
depending  on  the  potential  vector  and  potential  scalar  of  the 
y-subdomain,  respectively. 

The  inductive  term  of  (4)  is  given  by  the  expression: 

r 


f  /  jjnG(F,nds' 


subdomain  J,  respectively.  Taking  into  account  the  fact  that 
these  densities  are  constant  and  that  the  total  charge  in  the 
surfaces  is  equal  to  1/joi  we  have: 


=  f  V  -Lf  G{f,f')dl' 

47i:£ya)  “ 


—  f  G(r,r')dS'  dl 

o  M  J 


f  JJr')G(r,F')dS’  dl 


where  the  current  densities  on  the  junction  subdomain  j 
have  been  separated  into  two  parts  depending  on  which  part 
of  this  subdomain  we  are  considering:  J^{f')  on  the  surface 

of  the  wire  segment  (5j  )  and  7^(7')  on  the  surface  of  the 

BP  (5c“).  The  inductive  term  defined  in  (5)  can  be  largely 

J 

simplified:  a)  if  we  neglect  the  contribution  of  the  BP 
current  because  as  is  approximately  radial  its  equivalent 
momentum  vanishes;  b)  if  we  approximate  the  currentJ^ 
by  a  pulse  function  that  extends  over  the  first  half  of  the 
segment  of  the  junction  subdomain  j  and  c)  if  we  consider 
the  inductive  field  to  be  constant  along  the  wire  segment 
over  which  the  testing  function  of  the  i  subdomain  extends. 
The  simplified  expression  as  follows: 


Taking  into  account  that: 

b 

f  V/(r)dr  =f(b)-f(a} 

a 

expression  (8)  can  be  re-written  as: 


where  the  function  /  is  defined  as: 


f{l)  =  \  f  G{r,r)dl' 

h“  0 

-±fG(f.r,ds' 
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/  G{f^,r)di' 


where  is  the  centre  point  of  the  wire  segment 
corresponding  to  the  /-junction  subdomain.  On  the  other 
hand,  the  capacitive  term  of  (4)  can  be  obtained  as: 


Similar  expressions  for  the  rest  of  terms  of  the  MM  matrix 
can  be  obtained  considering  the  domain  and  the  limits  of 
integration  in  each  case.  Thus,  the  capacitive  and  inductive 

term  of  the  coefficient  J)  are  expressed  as: 
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/  V  /  p{f')G{f,f')dS' 

0  [ 

0  p** 
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ZZiiJ)  =  /  V  [w'jy)G{r,ndl'  dl  (12b) 


where  p(f')  and  o(f')  are  the  charge  densities  on  the 
wire  segment  and  on  the  BP  surface  of  the  junction 


where  Ij  is  the  wire  segment  of  the  y-antenna  subdomain 
and  J  represents  the  current  on  that  subdomain. 
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If  it  is  assumed  that  the  j  subdomain  of  the  body  is  shared 

by  the  Bezier’s  patches  A  and  B,  being  and  the 
current  of  the  basis  function  on  both  patches  respectively, 

the  coefficient  is  defined  as: 


-1 

4716/0) 


(13b) 


where  5^  and  are  the  areas  of  the  patches  A  and  B, 

Finally,  all  the  integral  involved  in  the  expressions  for  the 
coupling  impedances  are  carried  out  numerically  by  means 
of  the  Gauss’  Quadrature  Method.  This  algorithm  avoids  the 
singularity  of  the  Green’s  function  which  occurs  when  the 
points  f  and  f '  coincides.  Because  of  this,  the  expressions 
developed  above  hold  for  the  self-terms. 


from  [10],  The  sphere  has  been  modelled  by  a  mesh  of  BP 
that  form  4  meridians  and  9  parallel  lines  {76  basis 
functions  on  the  sphere).  Good  agreement  is  observed  in 
both  cases.  It  is  important  to  notice  that,  since  the  sphere  is 
a  body  of  revolution,  the  patches  on  the  poles  are  triangular 
and  considerably  smaller  than  the  remaining  patches  but  no 
special  treatment  of  the  current  has  been  applied. 


Figure  5.  Feeding  model  for  the  attached  antenna.  An 

impressed  voltage  of  0.5  Volts  is  imposed  on  the 
junction  subdomain  and  on  the  regular  surface 
subdomain  a,  b,  c  and  d  that  surrounds  the 
antenna 


In  the  study  presented  in  this  paper  the  magnetic  current 
frill  model  is  used  to  feed  the  antennas.  For  this  purpose,  a 
voltage  of  1  volt  is  distributed  among  the  subdomains  which 
share  the  BP  where  the  wire  is  attached,  as  shown  in  Fig. 5. 
This  figure  shows  the  digitized  impressed  voltage  which  is 
associated  to  each  one  of  these  subdomains  to  approximately 
model  the  magnetic  current  frill. 

3  RESULTS 

Several  geometries  have  been  analyzed  in  order  to  validate 
the  method  presented.  The  results  include  examples  of 
monopole  antennas  attached  to  perfectly  conducting  bodies 
such  as  spheres,  cylinders  and  other  ground  planes. 
Radiation  patterns  and  input  impedance  values  have  been 
obtained  and  compared  with  measurements  or  results  from 
other  methods. 

3.1  Monopole  on  conducting  sphere 


Theta  (degrees) 

—  Calculated  "  Calculated  [10] 

Figure  6.  Radiation  pattern  of  an  \/4  monopole  on  a  W4 
sphere 


Figs.6-7  shows  a  quarterwave  monopole  antenna  on  a 
conducting  sphere.  Computed  results  for  a^'K/4  and  a=W8, 
are  presented  (solid  line)  and  compared  with  those  taken 
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Monopole  on  an  polihedrical  conducting  body 


dB 


Figure  7.  Radiation  pattern  of  an  X/^  monopole  on  a 
sphere 

3.2  Monopole  on  a  conducting  cylinder 

Fig.  8  shows  the  configuration  of  a  monopole  antenna  (A) 
diametrically  opposite  a  passive  boom  (B)  on  the  surface  of 
a  cylinder  of  height  and  diameter  0.22m  and  0.20m 
respectively.  The  wavelength  is  0.36m.  Computed  values 
(solid  line)  for  the  <l)=0^  cut  are  compared  in  Fig.  8  with 
measured  data  (points)  taken  from  [8]  for  a  monopole  length 
=  0.08m  and  a  boom  length  =  0.44m.  Assuming  that  the  z- 
axis  coincides  with  the  cylinder  axis,  the  mesh  of  BP  of  the 
cylinder  is  defined  by  12  equidistant  meridian  lines  ((|)~ 
constant  curves)  and  5  parallel  lines  (z-constant  curves). 

dB 


Figure  8.  Radiation  pattern  of  a  monopole  antenna  (A) 
diametrically  opposite  a  passive  boom  (B) 


3.3 


A  very  special  conducting  box  is  shown  in  Fig.9.  A 
quarterwave  monopole  antenna  is  mounted  on  the  top  plate 
of  the  box  (the  most  narrow).  The  same  figure  also  includes 
a  sketch  of  the  mesh  considered  for  the  box  (eight  patches 
per  wavelength,  the  frequency  is  317  MHz).  The  total 
number  of  rooftops  on  the  body  is  1504.  Figs.  10-11 
illustrate  computed  (solid  line)  and  measured  (dashed  line) 
values  respectively  for  the  and  cuts,  respectively 
of  the  radiation  patterns.  It  is  important  to  notice  that, 
although  the  antenna  is  very  close  to  the  edges  of  the  model 
the  agreement  between  computed  and  measured  values  is 
relatively  good. 


1.10m 


(b) 

Figure  9.  (a)  Conducting  box  with  flat  sides;  (b)  Mesh  of 
BP  used  to  model  the  box 


dB 


Figure  10.  Radiation  pattern  of  a  monopole  antenna  on  the 
top  side  of  the  box  of  Fig.  9  (cut  <i>=(f) 
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dB 


—  Calculated  -  'Measured 

Figure  11.  Radiation  pattern  of  a  monopole  antenna  on  the 
top  side  of  the  box  of  Fig.  9  (cut 

3.4  Input  impedance  of  a  monopole  on  a  box 

Fig.  12  shows  a  monopole  antenna  centred  on  the  top  plate 
of  a  conducting  box.  The  computed  results  for  the  input 
impedance  of  this  antenna  are  compared  with  measurements 
taken  from  [11].  Several  densities  for  the  mesh  of  the 
conducting  box  have  been  considered  depending  on  the 
range  of  frequencies  considered  (a  number  of  ten  patches 
per  wavelength  is  always  guaranteed).  A  sketch  of  the  mesh 
for  the  frequency  band  1.75  to  4.0  GHz  is  included  in  the 
figure.  The  monopole  has  been  modelled  considering  20 
subdomains  for  all  the  frequency  band.  A  good  agreement 
between  calculated  and  measured  data  is  also  observed. 

3.5  Input  admittance  of  a  monopole  on  a  plane  with 
attached  reflected  plate 

Finally,  Fig.  13  shows  a  monopole  antenna  centred  on  a 
conducting  flat  plate  with  a  reflected  plate  attached  to  one 
side  of  the  ground  plane.  The  mesh  considered  for  all  the 
frequency  band  is  included  in  the  same  figure.  As  in  the 
previous  case  the  monopole  has  been  modelled  also  by  20 
subdomains.  The  radius  of  the  monopole  is  0.0008m.  The 
computed  results  for  the  input  admittance  of  this  antenna  are 
compared  with  measurements  taken  from  [7].  The  density 
for  the  mesh  of  the  conducting  plates  that  has  been 
considered  for  all  the  range  of  frequencies  analyzed  is 
shown  in  the  same  figure  (a  number  of  ten  patches  per 
wavelength  is  always  guaranteed).  Good  agreement  between 
computed  and  measured  data  is  observed. 


Frequency  (GHz) 


—  Measured  R  [11]  - '  Measured  X  [11]  •  Calculated  R  +  Calculated  X 

(b) 

Figure  12.  Monopole  antenna  on  a  conducting 

parallelepiped;  (a)  meshing  of  BP  for  the 
frequency  band  1. 75  to  4.0  GHz;  (b)  results  for 
the  input  impedance 

4  CONCLUSIONS 

A  new  technique  based  on  the  MM  to  analyze  antennas  on 
board  conducting  surfaces  has  been  presented.  The  scheme 
makes  use  of  modem  computational  geometry  tools  to 
model  the  scatterers  and  the  simplicity  of  wires  to  model  the 
on-board  antennas.  This  theory  has  been  implemented  in  a 
computer  code  in  order  to  obtain  radiation  patterns  or  input 
impedances  of  several  canonical  stmctures.  The  wire  can  be 
attached  to  any  patch  of  the  body  meshing  without  requiring 
any  special  attention.  Comparisons  between  numerical  data 
and  measurements  show  that  the  method  is  accurate  and 
efficient  for  radiation  patterns  analyses.  For  input  impedance 
the  approach  gives  reasonably  accurate  results. 
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Input  Admittance  (mmhos) 


—  Measured  G  [7]  “  ■  Measured  B  [7]  •  Calculated  G  +  Calculated  B 

(b) 


Figure  13.  Monopole  antenna  on  a  flat  ground  plane  with 
attached  reflected  plate;  (a)  meshing  of  BP  for 
all  the  frequency  range;  (b)  results  for  the  input 
admittance 
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Abstract 

A  brief  summary  of  the  variational  boundary- value  prob¬ 
lem  formulation  of  the  2D  finite  element/boundary  ele¬ 
ment  (FE/BE)  method  is  presented.  From  this  a  poste¬ 
riori  error  estimates  and  error  indicators  for  the  FE/BE 
method  are  developed  and  applied  to  electromagnetic 
scattering  and  radiation  problems.  The  results  obtained 
indicate  that  these  error  estimates  and  indicators  can 
be  obtained  within  negligible  computational  times  and 
can  be  used  successfully  to  obtain  valuable  a  posteriori 
accuracy  and  convergence  information  regarding  the  re¬ 
liability  of  the  FE/BE  method  solutions. 


1  Introduction 

The  2D  finite  element/boundary  element  (FE/BE) 
method^  has  been  used  extensively  over  the  past  few 
years  for  solving  electromagnetic  problems  numeri¬ 
cally  [1,  2,  3,  4,  5].  The  method  has  proved  to  be 
highly  successful  and  applicable  to  specifically  scattering 
and  radiation  problems  concerning  inhomogeneous,  ar¬ 
bitrarily  shaped  objects.  However,  the  solution  time  and 
memory  requirements  of  the  method  become  impractical 
when  electromagnetically  large  problems  are  considered. 
These  limitations  of  the  FE/BE  method  are  obviously 
dependent  on  the  computational  hardware  available. 

A  FE/BE  method  solution  results  in  approximated 
field  values  in  the  regions  under  consideration.  The  ac¬ 
curacy  of  these  solutions  is  dependent  on  the  specific 
problem  at  hand  as  well  as  the  approximation  functions 
used  with  the  FE/BE  method.  A  priori  knowledge  of 
the  accuracy  and  reliability  of  the  FE/BE  method  is  ob¬ 
viously  an  important  consideration  which  has  been  in¬ 
vestigated  by  a  number  of  researchers  [6,  1,7].  Another 

^The  technique  is  also  referred  to  as  the  finite  element /moment 
method  or  finite  element /integral  equation  method  in  the  litera¬ 
ture;  there  are  sometimes  differences  in  detail  but  the  basic  con¬ 
cepts  are  the  s£ime. 


important  consideration  would  be  an  a  posteriori  error 
estimate  of  FEM  or  FE/BE  method  solutions.  Reliable 
a  posteriori  error  estimates  could  serve  as  criteria  for  re¬ 
quired  accuracies  as  well  as  convergence  checks  for  the 
FE/BE  method  solutions  obtained.  A  posteriori  error 
estimates  for  the  FEM  as  well  as  the  BEM  have  been 
under  investigation  for  the  past  few  years  and  have  been 
applied  successfully  to  a  number  of  general  engineering 
problems  [8,  9,  10].  However,  very  little  has  been  pub¬ 
lished  on  this  topic  in  the  computational  electromagnetic 
literature  and  only  recently  has  work  started  appearing, 
for  example  [11,  12].  A  recent  special  issue  of  this  jour¬ 
nal  contained  several  papers  on  error  estimation,  none 
specifically  on  FEM  error  estimates  although  the  paper 
by  Hsiao  and  Kleinman  considers  boundary  integral  er¬ 
ror  control  [13].  Some  commercial  FEM  programs  in¬ 
clude  error  estimates  and  adaptive  meshes,  but  the  al¬ 
gorithms  are  often  proprietary. 

In  this  paper,  a  brief  summary  of  the  FE/BE  method 
formulation  will  be  presented  (section  2).  A  number  of  a 
posteriori  evTOT  estimates  and  error  indicators  for  FE/BE 
method  solutions  of  electromagnetic  problems  will  be 
formulated  in  section  3.  These  include  local  (in  each 
finite  element)  and  global  Element  Residual  Method 
(ERM)  error  estimates  [8,  9,  10],  T^-norm  boundary  field 
and  boundary  field  derivative  error  estimates  [14,  15],  a 
L^-norm  Neumann  boundary  condition  error  indicator 
and  a  radar  width  error  indicator.  It  will  be  shown  that 
these  are  highly  efficient  error  estimates  with  negligible 
computational  times  compared  to  the  solution  times  of 
FE/BE  method  solutions.  The  a  posteriori  error  esti¬ 
mates  and  error  indicators  developed  will  be  applied  to 
a  number  of  FE/BE  method  solutions  of  electromagnetic 
scattering  problems.  The  results  obtained  will  be  used 
to  investigate  the  accuracy,  reliability  and  applicability 
of  the  different  error  estimates  and  error  indicators  when 
applied  to  FE/BE  method  solutions  of  2D  electromag¬ 
netic  problems. 

In  section  4  general  conclusions  on  the  work  presented 
in  this  paper  will  be  drawn  and  further  research  that 
needs  to  be  done  will  be  discussed. 
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2  The  2D  FE/BE  Method 

In  this  section  the  2D  FE/BE  method  formulation  will 
be  presented^.  The  development  of  this  formulation  is 
not  new,  but  the  specific  notation  used  is  crucial  for  the 
development  of  the  FE/BE  method  error  estimates.  A 
iihorough  understanding  of  the  mathematics  behind  the 
FE/BE  method  is  essential  for  the  development  of  any 
kind  of  error  analysis  associated  with  the  method.  A 
rather  informal  discussion  of  the  functional  analysis  for¬ 
mulation  of  the  FE/BE  method  will  be  presented  in  this 
paper.  This  discussion  is  sufficient  to  prepare  the  reader 
for  the  development  of  the  error  estimates.  (Details  of 
the  formal  functional  analysis  formulation  of  the  FE/BE 
method  can  be  found  in  reference  [16].) 

The  formulation  will  be  presented  for  2D  Transverse 
Magnetic  (TM^)  polarized  scattering  problems.^  The 
TM^  scattering  formulation  is  a  special  case  of  the  more 
general  2D  scalar  Helmholtz  equation  formulation.  This 
equation  serves  as  the  governing  equation  for  a  variety 
of  2D  electromagnetic  problems,  including  static  elec¬ 
tric  and  magnetic  problems^,  as  well  as  closed  and  open 
boundary  electromagnetic  problems  for  TM^  and  TE^ 
polarization.  The  FE/BE  method  and  the  a  posteri¬ 
ori  error  estimate  methods  considered  in  this  paper  are 
applicable  to  all  the  above  mentioned  electromagnetic 
problems. 


2,1  Variational  boundary- value  problem 
formulation  of  the  FEM 

Consider  figure  1.  For  2D  scattering  problems  involving 
a  TM^  polarization  incident  field,  the  2D  electromag¬ 
netic  field  equation  can  be  written  as  a  boundary  value 
problem  [6,  pp.l85]  [1,  pp. 72-73].  (Throughout  this  pa¬ 
per  the  electric  field  component  associated  with  the  TM^ 
polarized  field,  Ez^  will  be  denoted  by  E): 

V—VE  +  €rklE  =  0  in  ft  (1) 

f^r 

^  Only  two-dimensional  problems  will  be  considered  in  this  pa¬ 
per.  Variation  in  the  2D  x  —  y  plane  will  be  assumed  (for  all  fields 
and  structures)  with  no  variation  in  the  ^-direction  (extending  to 
infinity). 

^The  convention  for  polarization  in  a  2D  plane  which  will  be 
used  throughout  this  paper  is  as  follows:  For  TM^  polarization 
the  magnetic  field  vector  has  only  x  and  y  components,  thus  f} 
is  always  transverse  to  the  z  pleme.  In  this  case  H  has  only  a  z 
component.  For  TEz  polarization  the  electric  field  vector  has  only 
X  and  y  components,  thus  S  is  ^dways  transverse  to  the  z  plane. 
In  this  case  A  has  only  a  z  component. 

*In  these  cases,  the  2D  scalar  Hehnholtz  equation  reduces  to 
the  2D  Laplace  or  Poisons  equations. 


with  Dirichlet  boundary  condition  on  a  perfectly  con¬ 
ducting  region 

E  =  0  on 

(2) 

and  Neumann  boundary  condition  on  the  fictitious 
boundary 

1^=51  on  (3) 

on 

where  gi  are  the  prescribed  values  for  on  .  The 
boundary- value  problem  of  equations  1  to  3  can  be  writ¬ 
ten  as  a  variational  bound  ary- value  problem  [6,  pp.219]: 

/  (VE  •  —  Vv  -  erklEv)  dft  =  -  /  giv  dV'^^  (4) 
Jn  Jr^i 

with  V  an  arbitrary  weighting  or  testing  function  on 
n.  Using  conventional  finite  element  procedures  [17][16, 
pp.l4],  equation  4  can  be  written  in  matrix  form: 

[S][E]  +  [T][^]  =  0  (5) 

with  [5]  the  FEM  system  matrix,  [E]  the  unknown  field 
coefficient  column  matrix,  [T]  the  FEM  boundary  condi¬ 
tion  system  matrix  and  [|^]  the  field  derivative  column 
matrix.  Equation  5  is  the  FEM  matrix  equation  for  TM^ 
polarized,  2D  electromagnetic  problems.  The  solution 
of  the  FEM  matrix  equation  5  yields  the  approximated 
field  solution,  E^  if  gi  of  equation  3  is  known  and  used 
to  construct  the  column  matrix  [§“]• 


Figure  1:  Finite  region,  fl,  in  which  the  FEM  will  be 
applied,  enclosed  by  a  boundary,  .  A  is  an  exterior 
free-space  region  extending  to  infinity. 


2.2  The  boundary  element  method 

Many  practical  electromagnetic  field  problems,  such  as 
scattering  and  radiation  problems,  are  open  boundary 
problems.  The  domain  of  these  problems  extends  to  in¬ 
finity,  making  the  finite  element  method  unsuitable  for 
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the  solution  of  such  problems.  One  way  to  overcome  this 
problem  is  to  couple  the  finite  element  method  with  an 
unbounded  region  method  such  as  the  boundary  element 
method.  In  this  way,  the  advantages  of  both  the  FEM 
and  BEM  are  exploited  [1][5]. 

A  fictitious  boundary  denoted  by  can  be  defined 
around  the  structure  under  consideration  dividing  the 
solution  region  into  an  interior  region,  fi,  and  an  exterior 
region,  A  (see  figure  1).  A  BEM  solution  of  the  normal 
field  derivatives  on  can  be  used  as  Neumann  bound¬ 
ary  condition  for  the  boundary  value  field  problem  in  Q. 
The  FEM  can  thus  be  applied  to  the  field  problem  in 
n  with  appropriate  boundary  conditions  (obtained  with 
the  BEM  solution)  on  . 

Consider  the  exterior  region  A  bounded  by  F^*  and 
infinity.  For  electromagnetic  field  problems,  the  homo¬ 
geneous  free-space  Helmholtz  equation  must  be  satis¬ 
fied  in  this  region.  The  variational  form  of  the  bound¬ 
ary  element  method  integral  equation  [1,  pp.299][14, 
pp.41][16,  pp.l7]  can  be  written  in  the  BEM  matrix 
equation  form  [16,  pp.l9]: 


[H][E]-[G\ 


dE 

dn' 


(6) 


The  matrix  elements  of  [H]y  [G]  and  can  be  calcu¬ 

lated  using  the  boundary  element  integral  equation  [16, 
pp.l7]  and  the  matrix  elements  of  [E]  and  [|^]  ^.re  un¬ 
known  field  coefficient  values  which  can  be  obtained  with 
the  BEM  solution.  The  elements  of  matrix  [F*”^]  contain 
external  source  information.  In  the  rest  of  this  paper,  E 
represents  the  FEM  and  E  the  BEM  computed  solution. 


2.3  Coupling  the  FEM  and  BEM 


For  a  closed  region,  Q,  a  FEM  solution  is  obtain¬ 
able  if  certain  boundary  conditions  are  known.  These 
could  be  Dirichlet  and/or  Neumann  boundary  condi¬ 
tions.  For  open  boundary  problems  one  can  use  the 
fictitious  boundary  (discussed  in  the  previous  sections) 
as  closure  of  the  finite,  interior  region,  Q,  with  inhomo¬ 
geneous  Neumann  boundary  conditions  specified  on  the 
fictitious  boundary  .  By  using  ^  of  equation  6  as 
the  Neumann  boundary  condition  on  F^^ ,  a  solution  of 
the  FEM  in  Q  is  obtainable.  That  is,  gi  of  equation  3 
can  be  set  equal  to  of  equation  6  on  F^^ .  (The  mi¬ 
nus  sign  is  due  to  the  direction  of  the  normal  n  —  see 
figure  1): 


dE 


dn' 


on 


(7) 


Using  the  above  equation  together  with  equations  4 
and  5  leads  to  a  modified  FEM  matrix  equation 

[5][^]  -  [T] 

Setting  E  of  equation  5  equal  to  E  of  equation  6  on  F^^ 
and  coupling  the  two  matrix  equations  6  and  8  (each  con¬ 
taining  two  unknown  coefficient  column  matrices)  is  thus 
equivalent  to  solving  the  FEM  matrix  equation  of  equa¬ 
tion  5  with  the  Neumann  boundary  condition  of  equa¬ 
tion  7  on  F^^ . 


dE 


dn* 


==0 


(8) 


3  A  Posterior  Error  Estimates 

A  reliable  a  ’posteriori  error  estimate  for  a  FE/BE 
method  solution  would  enable  one  to  obtain  valuable 
convergence  information  without  having  to  solve  the 
same  problem,  using  a  larger  number  of  unknowns.  The 
error  estimate  can  thus  be  used  as  a  convergence  check 
for  practical  electromagnetic  problems  for  which  no  an¬ 
alytical  solution  exists.  This  would  be  especially  advan¬ 
tageous  when  electromagnetically  large  problems,  reach¬ 
ing  the  practical  solution  time  and  memory  limits  of  the 
computer  at  hand,  are  considered.  This  is  true  regardless 
of  the  computational  power  available. 

Local  and  global  finite  element  method  error  estimates 
have  been  under  investigation  for  the  past  decade  [8, 
9,  10].  Of  these,  the  Element  Residual  error  estimate 
Method  (ERM)  seems  to  be  the  most  reliable  and  the 
easiest  implementable  method.  Applications  of  this 
FEM  error  estimate  method  were  directed  mainly  at 
mechanical  and  fluid-dynamical  problems  although  the 
method  is  applicable  to  a  wide  variety  of  positive  defi¬ 
nite  FEM  type  problems  [18].  The  governing  equation 
for  static  electric  field  problems  (Laplace’s  or  Poisson’s 
equation)  is  positive  definite,  and  the  successful  appli¬ 
cation  of  the  ERM  to  such  problems  has  recently  been 
published  [19].  The  ERM  applied  to  wave  equation  type 
problems,  of  which  electromagnetic  scattering  and  radi¬ 
ation  problems  are  examples,  will  be  considered  in  this 
section.  Although  the  FEM  system  equations  for  such 
problems  are  not  necessarily  positive  definite,  an  element 
residual  type  error  estimate  method  will  be  introduced, 
which  is  applicable  to  such  problems. 

A  L^-norm  boundary  field  error  estimate,  which  has 
been  applied  to  linear  acoustic  problems  [14,  15],  and  a 
L^-norm  boundary  field  derivative  error  estimate  have 
been  investigated  and  will  be  modified  for  application  to 
the  BEM  part  of  the  FE/BE  method  solution.  A  L^- 
norm  Neumann  boundary  condition  error  indicator 
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for  the  FE/BE  method  will  also  be  introduced.  The 
norm  error  estimates  and  error  indicator  will  be  used 
to  develop  far-field  error  indicators  for  FE/BE  method 
solutions  of  electromagnetic  problems.  These  far-held 
error  indicators  are  especially  important  when  consid¬ 
ering  electromagnetic  problems  concerning  radar-width 
and  radiation  intensity  results,  obtained  with  the  FE/BE 
method  solution. 

The  mathematical  development  of  the  error  estimates 
and  error  indicators  presented  in  this  section,  is  based 
on  the  functional  analysis  formulation  of  the  FE/BE 
method,  presented  in  section  2  and,  in  more  detail,  in 
reference  [16]. 

3,1  Error  Estimates  for  the  FEM 

The  Element  Residual  error  estimate  Method  (ERM)  for 
FE/BE  method  solutions  of  electromagnetic  scattering 
problems  will  be  formulated  in  this  section.  This  is  an 
a  posteriori  error  estimate  method,  and  can  be  applied 
after  a  FE/BE  method  solution  has  been  obtained.  Let 
us  denote  a  “hrst”  approximated  electric  held  solution 
obtained  in  Q  as  and  a  “second”  (higher  order)  held 
solution  as  E^.  (Note  that  fi  is  the  closure  of  Q,  i.e. 
the  domain  Q  and  its  boundary  F^').  The  ERM  will  be 
applied  after  solution  E^  has  been  calculated  and  yields 
an  estimate  of  the  relative  FE/BEM  electric  held  error: 

£-12  ^  £.2  _  £>1  (9) 

The  above  mentioned  error  will  be  quantihed  in  terms 
of  a  norm  (as  required  by  the  ERM).  For  static  elec¬ 
tric  held  problems  an  appropriate  norm  is  the  energy 
norm  [19]  which  yields  a  quantitative  error  value  of  the 
stored  energy  in  the  region  under  consideration.  For 
electromagnetic  problems  a  similar  norm  associated  [16, 
pp.248]  with  the  stored  electric  and  magnetic  energy  in 
the  closed  region  Q  will  be  employed.  This  norm  will  be 
called  the  electromagnetic  energy  norm  (EM-norm)  and 
is  given  by  [16,  pp.249]: 

M 

ifc  =  l 

with  M  the  number  of  hnite  elements  and 

Jctk  \ 

Re{erko}\E^^\‘^)  dQk  (H) 
the  local  EM-norm  associated  with  finite  element  Clk- 

The  ERM  will  thus  yield  a  quantitative  error  estimate 
of  the  stored  electric  and  magnetic  energy  in  the  region 


Q.  This  error  estimate  will  be  denoted  by  • 

The  EM-norm  satishes  the  requirements  for  a  norm  gen¬ 
erated  by  an  inner  product  [16,  pp.246]  and  is  not 
unique  for  ERM  error  estimates  applied  to  electromag¬ 
netic  problems.  Other  norms  could  be  developed  and 
might  lead  to  improved  error  estimate  results.  It  will 
however  be  shown  that  the  EM-norm  introduced  here  re¬ 
sults  in  reliable  ERM  error  estimates  for  FE/BE  method 
solutions  of  general  (2D)  electromagnetic  problems  (in¬ 
cluding  problems  involving  lossy  materials). 


3.1.1  Local  ERM  error  estimates 


An  estimate,  of  actual  error 

in  each  finite  element,  Qk,  can  be  ob¬ 
tained  by  solving  the  following  local  problem  [10, 
pp.l29][20][16],  related  to  the  EM-norm  (a  sub-  or  su¬ 
perscript  k  denotes  that  a  quantity  is  associated  with 
finite  element  Qa;): 


/  d^lk  = 

I 


dji 

duk 


(TCk  for  all 


(12) 


In  this  equation  vl^  is  the  weighting  functions  on  Q. 
present  in  the  FE/BEM  solution  but  not  in  the  solu¬ 
tion  E^.  Tk(i)  is  the  side  of  the  triangular  finite  element 
Qk  connecting  Qk  with  its  /*'•  neighboring  element 
Also,  rl  is  the  local  element  residual  of  the  governing 
equation  (equation  1)  cissociated  with  E^,  with 


4 


1 


VEl-elklEl 


(13) 


is  the  jump  in  normal  derivative  between  el¬ 
ement  Q.k  and  its  neighboring  element,  element  Q,k{i), 
with 

f  dEj  \  ^  (14) 

I  drikQ)  J  dnk(i)  /r?  dnk(j)  j 


Equation  12  appears  formidable  but  can  be  written  in 
the  following  matrix  form  for  each  finite  element,  Q*  [16, 
pp.l03]; 


_  JJ7’*:(12)j  (15) 

with  a  square  matrix  with  matrix  elements  which 

can  be  calculated  using  the  integral  on  the  left  hand  side 
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of  equation  12;  is  a  column  matrix  with  ma¬ 

trix  elements  the  unknown  local  error  estimate  coeffi¬ 
cients;  is  a  column  matrix  with  matrix  elements 

which  can  be  calculated  using  the  integrals  on  the  right 
hand  side  of  equation  12.  The  elements  of  can 

be  computed  rapidly  using  tables  for  the  result  of  in¬ 
tegration  over  a  prototypical  element  along  the  lines  of 
[17,  Chapter  4].  (Details  are  given  in  [16,  Appendix  F]). 
Solving  the  local  matrix  equation,  equation  15,  is  equiv¬ 
alent  to  solving  the  local  element  residual  error  estimate 
equation,  equation  12,  yielding  and  thus  the  local 
error  estimate  Note  that  equation  15  is  a 

very  “small”  matrix  equation  involving  only  one  finite 
element  (at  a  time).  This  equation  can  thus  be  solved 
within  negligible  computational  time. 

The  specific  form  of  equation  12  is  important  to  ensure 
the  following  relation  between  the  global  error  estimate, 
and  actual  global  error  in  the  EM-norm  [20] [10, 
pp.l26,ppl28]: 

l|^12||gM  ^  ^12||^12||gM  Qg) 

with  the  global  error  estimate  inequality  constant 
for  the  error  estimate  in  the  EM-norm.  It  is  obvious 
that  the  accuracy  of  the  global  error  estimate  depends  on 
the  value  of  the  constant  in  equation  16.  is  de¬ 
pendent  on  the  order  of  the  basis  functions  used  as  well 
as  the  shape  of  the  finite  elements.  In  practice  =  1 
is  assumed  with  all  error  estimates.  A  better  approxi¬ 
mation  of  would,  however,  improve  the  accuracy  of 
the  error  estimates  [10,  pp.l29][8]. 

A  discussion  of  the  contribution  of  each  term  of  equa¬ 
tion  12  to  the  error  estimate  will  now  follow.  This  should 
give  a  qualitative  understanding  of  the  ERM  and  specif¬ 
ically  its  application  to  EM  problems. 

The  weighted  residual  integral  term 

The  first  term  on  the  right  hand  side  of  equation  12 
is  a  weighted  residual  term.  The  “first”  solution  re¬ 
sults  in  a  field  solution  El  in  element  Qk-  The  residual 
rl  will  be  zero  everywhere  in  fijt,  if  E^  is  the  true  (or 
exact)  field  solution.  This  is  due  to  the  fact  that  equa¬ 
tion  13,  with  rl  =  0,  is  a  form  of  the  governing  equation 
which  will  be  satisfied  by  the  true  field  solution,  E^ 
is  an  approximate  finite  element  solution  the  residual  rl 
will  be  non-zero  and  can,  in  general,  vary  in  numerical 
value  over  the  finite  element  .  The  numerical  value  of 
the  weighting  functions,  up,  on  Qk  is  zero  at  all  nodal 
points  on  Qk  associated  with  the  basis  functions  of  the 
solution  E^.  The  numerical  value  of  on  fljb  is  one  at 
all  the  nodal  points  on  Qk  associated  with  the  solution 
E^,  but  not  with  the  solution  in  E^ .  At  non-nodal  points 
vl^  is  a  value  between  zero  and  one  varying  to  the  order 


of  the  basis  functions  used  with  E^^  Multiplication  of 
the  residual  rl  with  vl^  and  integration  over  is  thus 
equivalent  to  weighting  the  residual  with  over  Qk 
This  weighted  term  will  thus  have  a  relatively  large  nu¬ 
merical  value  if  the  residual  is,  on  average,  large  in  the 
regions  where  the  “new”  nodal  points,  corresponding  to 
the  “new”  basis  functions  of  E^  on  fijb,  will  be  situated. 
The  weighted  term  will  have  a  small  numerical  value  if 
the  residual  is,  on  average,  small  in  the  regions  where 
the  “new”  nodal  points,  corresponding  to  the  “new”  ba¬ 
sis  functions  of  E^  on  Q*,  will  be  situated. 

Weighted  normal  derivative  continuity  integral 
term 

The  second  integral  term  on  the  right  hand  side  of 
equation  12  is  a  weighted  normal  derivative  continuity 
term  on  '^k{i)-  If  is  the  true  field  solution,  continu- 

1  dE^ 

ity  between  the  terms  and  of  equa¬ 

tion  14  must  hold  on  '^k(i)-  This  continuity  requirement 
is  related  to  the  zero  divergence  requirement  [21,  pp.l27]. 
This  means  that  {9^}  of  equation  14  must  be  zero 

for  a  true  field  solution.  The  different  local  basis 
functions  (or  approximation  functions)  in  each  finite  el¬ 
ement  Qk  result  in  discontinuities  on  'Tk(i)  for  E^  an 
approximate  finite  element  solution.  In  the  integral  of 
equation  12,  general,  vary  on  Tfc(/).  The 

term  |  \  is  weighted  by  the  restriction  of  the  vP 

to  Tfc(/). 

Weighted  Neumann  boundary  condition  inte¬ 
gral  term 

The  third  integral  on  the  right  hand  side  of  equation  12 
is  a  weighted  Neumann  boundary  condition  integral  on 
Tfc  .  The  Neumann  boundary  condition  will  be  sat¬ 
isfied  exactly  with  E^  the  true  field  solution,  making  the 
quantity  (^gi  -  zero.  If  E^  is  an  approx¬ 

imate  finite  element  solution,  the  Neumann  boundary 
condition  on  Tjb  H  will  be  satisfied  approximately 

and  the  term  (gi-  )  will  be  non-zero.  In 

iTfcor^i 

the  integral  of  equation  12,  (31  - 

general,  vary  on  .  The  term  ^5-1  - 

is  weighted  with  the  restriction  of  to  Tjb  O  . 

3.2  Error  Estimates  for  the  BEM 

In  this  section,  two  L^-norm  [6,  pp.50,72]  error  estimates 
for  the  FE/BE  method  solution,  E^ ,  will  be  introduced. 
Both  error  estimates  are  associated  with  the  accuracy 
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of  the  field  solution  on  ,  and  are  related  to  global 
field  and  field  derivative  errors  on  .  These  error  es¬ 
timates  are  estimates  of  the  ‘‘true”  errors;  that  is  the 
errors  in  relative  to  the  true  field  solution  [14,  15]. 
Boundary  errors  appear  to  have  been  seldom  used  in 
computational  electromagnetics  as  a  measure  of  solution 
accuracy;  the  only  example  we  are  aware  of  is  [22],  which 
includes  plots  of  the  tangential  electric  (error)  field  along 
a  wire  antenna. 


3,2.1  The  boundary  field  L^-norm  residual  error 
estimate 


The  first  L^-norm  error  estimate  is  associated  with  the 
accuracy  with  which  the  boundary  element  approxima¬ 
tion  functions  can  approximate  the  actual  fields  on  the 
boundary  The  FE/BE  method  solution  yields  field 
and  field  derivative  coefficients  which,  together  with  the 
boundary  element  basis  functions,  can  be  used  to  ap¬ 
proximate  the  fields  and  field  derivatives  on  .  These 
field  approximations  can  be  used  to  calculate  the  fields 
at  any  point  in  the  exterior  region  A,  using  a  numerical 
approximation  of  Huygens’  principle  [16,  pp.57].  The 
field  approximations  can  also  be  used  to  calculate  the 
fields  on  F^^  using  the  following  equation  [14,  pp.40]; 


^1 


d^{fo,rs) 

dn'(fs) 


dn'(fs) 


'if(ro,rs) 


dr'^'ifs) 


+2F‘'^%ro) 


(17) 


On  the  boundary,  the  FE/BE  formulation  ensures  that 
E^(fs)  =  E^(fs)  and  they  can  be  used  interchangeably 
here.  The  quantities  To  fg  are  the  BEM  observation 
and  source  points  on  F^^,  is  2D  homoge¬ 

neous  Green’s  function  and  E^^^(ro)  is  the  incident  elec¬ 
tromagnetic  field  value  at  fo .  E^(ro)  is  the  approximated 
field  on  calculated  using  equation  17. 

A  boundary  field  residual  can  now  be  defined  as^: 

Ri.^,(fo)  =  FHro)-F'(ro)  (18) 


with  iJpN,  (n,)  the  BEM  residual  for  the  FE/BE  method 
field  solution,  F^(fc),  at  fg,  on 

^Note  that  E^(fo)  is  the  approximate  field  solution  (at  Tq  on 
obtained  with  the  “first”  FE/BE  solution.  The  function 
(fo)  is  thus  of  order  equal  to  the  order  of  the  FEM  and  BEM 
basis  functions.  The  function  E^  (fo)  is  also  an  approximated  field 
solution  (at  fo  on  F^i)  obtained  with  equation  17  (and  thus  in¬ 
directly  from  the  “first”  FE/BE  method  solution).  The  function 
E^  (fo)  is  of  much  higher  order  than  E^  (f©)  due  to  the  presence 
of  the  Green’s  function,  ^'(fo,fs)j  ^nd  the  incident  field  function, 
E^^^(ro),  in  equation  17. 


A  global  boundary  field  L^-norm  residual  error  esti¬ 
mate  can  be  obtained  using  equation  18: 

\\e],.A\h  =  (19) 

It  is  evident  that  the  residual  (r^)  will  be  zero  if 
is  the  true  field  solution  and  is  the  true  field 

derivative  solution  on  for  then  the  approximate 

fields,  E^{fo)  and  E^{fo)  will  both  be  equal  to  the  true 
field  solution  at  r^.  From  equation  19  it  is  also  evident 
that  the  boundary  field  error  estimate,  HGp^i  11^2,  will  be 
zero  if  E^  is  the  true  solution.  Asymptotic  exactness  and 
upper  and  lower  bounds  for  the  error  estimate  l|0pNi  |1^2 
have  not  yet  been  obtained  for  general  FE/BE  method 
solutions  (although  this  has  been  done  for  some  special 
cases  [14,  pp.56])  mainly  due  to  the  difficulties  arising 
from  the  non-local  nature  of  the  BEM  when  applied  to 
general  electromagnetic  problems. 

Typically,  this  error  estimate,  which  indicates  the  de¬ 
gree  of  accuracy  of  the  BEM  part  of  the  solution,  con¬ 
verges  more  rapidly  than  the  FEM  error  estimates,  as 
will  be  seen  in  the  examples  given  in  §3.4.  However,  it 
remains  a  useful  estimate  and  is  required  for  the  far-field 
error  indicators  to  be  developed  in  §3.3. 


3.2.2  The  boundary  field  derivative  L^-norm 
residual  error  estimates 

Procedures  similar  to  those  described  in  the  previous 
subsection  can  be  followed  to  obtain  a  BEM  residual 
for  the  FE/BE  method  field  derivative  solution  [16, 

pp.106]: 

61  dE\fg) 

(20) 


is  the  normal  derivative  of  E^{fo)  of  equa- 

dn  (fo)  •  '  1 

tion  17.  The  first  term  on  the  right  of  equation  20  is  the 
field  derivative  computed  using  the  BEM  part  of  the  so¬ 
lution  (and  has  the  same  order  of  accuracy  as  the  FEM 
part,  since  the  BEM  expands  both  the  field  and  field 
derivative  to  the  same  order);  the  second  term  is  com¬ 
puted  numerically  using  the  Green’s  function  of  equa¬ 
tion  17. 

A  residual  for  the  FE/BE  method  Neumann  bound¬ 
ary  condition  is  [16,  pp.108]: 


(^o)  — 


dE\fo) 

dn{ro) 


dE^rp) 

dn{fo) 


(21) 
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This  estimate  differs  subtly  from  that  in  equation  20, 
in  that  the  residual  -RpAr,  (r^)  is  the  difference  (at  fo,  and 
on  r^‘)  between  the  normal  field  derivative  of  £'^(fo), 
obtained  by  numerically  differentiating  the  FEM  part  of 
the  solution,  and  the  normal  field  derivative  obtained  ex¬ 
plicitly  from  the  BEM  part  of  the  solution  (whereas  the 
estimate  of  equation  20  is  the  difference  between  the  so¬ 
lution  obtained  directly  with  the  BEM  and  that  obtained 
indirectly  from  the  BEM  via  the  Green’s  function  inte¬ 
gral  of  equation  17).  The  terms  in  equation  21  differ  due 
to  the  indirect  enforcement  of  the  Neumann  boundary 
conditions  —  see  section  2.3  —  and  the  different  orders 
of  approximation  functions  (the  field  derivative  solution 
computed  with  the  FEM  is  of  one  order  less  than  that 
computed  with  the  BEM). 

Explicit  numerical  results  for  these  estimates  will  not 
be  given  in  this  paper,  but  these  estimates  are  used  in 
§3.3  and  implicitly  contained  in  the  results  presented  in 
§3.4  and  have  thus  been  defined  here. 


3.3  Far-field  error  indicators 


The  FE/BE  method  solution  results  in  approximate  val¬ 
ues  for  the  field  and  field  derivatives  on  the  boundary 
.  From  this  one  can  calculate  the  field  value  at  any 
point  (in  the  near,  intermediate  or  far-field)  in  the  exte¬ 
rior  region.  A,  using  a  numerical  approximation  of  Huy¬ 
gens’  principle  [16,  pp.57]  (the  exterior  BEM  equation). 

A  far-field  residual  quantity  can  be  defined  as; 


R^(ro) 


(^5) 


g^(ro,f,) 

dn'(fs) 


(22) 


where  the  residuals  of  equations  18  and  21  have  been 
substituted  into  the  BEM  equation  in  the  place  of 

(r,)  of  equation  20  can  be 

used  instead  of  Rpui(^s),  but  the  indicator  would  then 
be  related  to  the  accuracy  of  only  the  BEM  part  of  the 
FE/BE  method  solution. 

The  far-field  residual  R^(fo)  will  be  zero  if  is  the 
true  solution,  for  then  both  residuals,  .RpA^j(r5)  and 
(^a),  used  in  the  above  equation,  will  be  zero. 

A  radar-width  error  indicator  can  now  be  calculated 
using  the  far-field  residual,  R^(fo)^  in  the  radar-width 
equation  [16,  pp.57]  instead  of 


y((/))  =  27r|fo| 


J^inc 


(23) 


A  radiation  intensity  error  indicator  can  also  be  calcu¬ 
lated  using  the  far-field  residual,  R^{fo),  in  the  radiation 
intensity  equation  [16,  pp.l09]  instead  of 

k{<l>)  =  2'K\fo\\R}{ro)^  (24) 

These  error  indicators  estimate  the  possible  error  in  the 
radar-width  and  radiation  intensity  due  to  the  errors 
in  the  approximate  FE/BE  method  field,  and  field 
derivative  solution.  Both  these  error  indicators  will  be 
zero  if  is  the  true  field  solution. 


3.4  Numerical  Results 

The  error  estimates  and  error  indicators  formulated  in 
the  previous  sections  have  been  applied  to  a  variety  of 
electromagnetic  problems,  with  promising  success  [16]. 
Two  examples  will  be  considered  in  this  section.  It  will 
be  shown  that  a  combination  of  the  different  error  esti¬ 
mates  and  indicators  can  be  used  as  a  reliable  a  posteri¬ 
ori  estimate  of  the  accuracy  and  convergence  of  FE/BE 
method  solutions. 


3.4.1  A  Scattering  Example 

The  first  numerical  example  which  will  be  considered  is 
for  electromagnetic  scattering  from  the  lossy  (using  prac¬ 
tical  material  parameters)  dielectric  right-circular  cylin¬ 
der  of  figure  2  excited  by  a  TE^  polarized,  plane  incident 
field  at  a  frequency  of  2  GHz.  (The  TE^  error  estimates 
are  analogous  to  the  TM^  estimates  given  in  the  preced¬ 
ing  sections).  Local  error  estimates  and  actual  errors  are 
compared  with  each  other  in  figures  4  and  5®.  In  table  1 
percentage  global  estimated  errors  (in  the  appropriate 
norms)  are  compared  to  actual  percentage  errors.  Note 
that  in  tables  1  and  3,  percentage  estimated  and  actual 
errors  of  more  than  100%  are  denoted  by  >100%.  Per¬ 
centage  error  estimates  are,  in  general,  not  much  larger 
than  100%,  even  if  the  actual  percentage  error  is  much 
larger  than  100%.  A  percentage  error  near  100%  does, 
however,  indicate  that  the  solution  at  hand  is  not  at 
all  close  to  the  true,  or  converged,  solution.  It  is  evi¬ 
dent  that  the  global  EM-norm  error  estimate  is  a  good 
estimate  of  the  actual  error.  The  global  percentage  er¬ 
ror  estimates  compare  reasonably  well  with  the  actual 
global  percentage  errors.  Important  to  note  is  that  the 

®The  finite  element  numbers  in  all  the  “local  EM-norm  error 
comparison”  figures  in  this  section  have  been  arranged  in  such 
a  way  that  the  corresponding  actual  EM-norm  error  values  are 
presented  in  increasing  order.  The  error  estimate  values  plotted 
are  local  EM-norm  error  estimate  values  in  the  corresponding  finite 
elements. 
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global  EM-norm  error  estimate  clearly  identifies  the  so¬ 
lutions  which  are  not  at  all  close  to  convergence,  and 
gives  a  reasonably  accurate  estimate  of  the  percentage 
errors  for  the  solutions  which  are  close  to  convergence  or 
have  converged  to  the  true  solution. 


Figure  2:  Double  dielectric  right-circular  cylinder 
(tri  =  1.2  + jl.5,  er2  =  1.5  + j2,  ri  =  0.27n,  r2  =  0.3m). 
The  cylinder  is  excited  by  a  plane  wave  incident  field 
with  TEz  polarization.  (Note  that  the  time  conven¬ 
tion  is 

A  “dB”  error  margin  can  defined  as: 

=  10\log{'f{<l))}  -  log{j{(f))  -  f(</’)}l  (25) 

where  is  the  radar- width  calculated  using  the 

FE/BE  method  and  is  the  far-held  error  indica¬ 
tor  of  equation  23.  Notice  that  and  thus  Aj{<j)) 

will  be  zero  if  is  the  true  field  solution. 

The  “dB”  error  margin  for  the  first  and  second  or¬ 
der  FE/BE  method  solutions  of  the  radar- width  of  the 
scattering  problem  under  consideration  are  shown  in  fig¬ 
ure  6.  The  first  and  second  FE/BE  method  solution  as 
well  as  the  analytical  solution  of  the  radar-width  are  also 
shown. 

All  error  estimates  and  error  indicators  indicate  that 
the  second  order  basis  function  FE/BE  method  solution 
for  M  =  1502  (M/A2  =  67)  and  =  90  (Mj/A  5.4) 
has  converged  to  an  acceptably  accurate  solution  (a  max¬ 
imum  field  error  in  the  FEM  region  of  around  5%  and  in 
the  radar-width  of  O.ldB).  This  is  confirmed  by  the  ac¬ 
tual  errors  and  the  radar-width  results  presented.  These 
error  estimates  and  error  indicators  have  been  obtained 
within  negligible  computational  times  (see  table  2)  com¬ 
pared  to  the  computationally  expensive  FE/BE  method 
solutions.  This  is  especially  true  for  the  FE/BE  method 


Figure  3:  Radiating  2D  horn  antenna,  with  a  =  0.55m, 
b  =  0.2m,  c  =  0.45m,  d  =  0.35m,  e  =  0.2m  and  /  =  0.05m. 


solutions  concerning  second  and  third  order  basis  func¬ 
tions  and  relatively  large  numbers  of  finite  and  boundary 
elements. 

3.4.2  A  Radiation  Example 

The  second  numerical  example  which  will  be  considered 
is  for  electromagnetic  radiation  from  the  2D  horn  an¬ 
tenna  of  figure  3.  The  horn  antenna  is  excited  by  the 
TEi  mode  aperture  field  (frequency:  750  MHz): 

E,  =  Eocos{^)  (26) 

Local  error  estimates  and  actual  errors  are  compared 
with  each  other  in  figures  7  and  8  and  it  is  again  evi¬ 
dent  that  these  local  error  estimates  clearly  identify  the 
regions  where  the  largest  EM-norm  errors  occur.  In  ta¬ 
ble  3,  percentage  estimated  errors  (in  the  appropriate 
norms)  are  compared  to  actual  percentage  errors.  The 
L^-norm  boundary  field  error  estimates  indicate  that 
enough  boundary  elements  are  used  to  approximate  the 
fields  on  the  boundary  for  all  values  of  considered. 
This  is  even  true  for  =  35  (M^/A  =  4.5)  with  first 
order  boundary  element  basis  functions.  (Note  that  in 
this  case  the  error  estimates  show  that  enough  bound¬ 
ary  elements  were  used.  The  actual  errors  are  quite  large 
due  to  the  fact  that  too  few  finite  elements  were  used). 
The  EM-norm  ERM  error  estimates  indicate  that  the 
second  order  basis  function  FE/BE  method  solution  for 
M  =  859  (M/A  =  145)  and  Mb  =  76  (M^/A  =  9.8)  has 
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converged  to  a  satisfactorily  accurate  FE/BE  method 
solution.  (This  corresponds  to  around  100  elements  per 
square  wavelength).  This  is  confirmed  by  the  actual  EM- 
norm  errors. 

A  MB”  error  margin  for  the  FE/BE  method  solution 
of  the  radiation  intensity  can  be  obtained  using  the  radi¬ 
ation  intensity  error  indicator  of  equation  24.  This  error 
margin,  denoted  by  AA"(0),  can  be  calculated  as  (this  is 
similar  to  the  radar- width  error  margin  of  equation  25); 

^k{<j>)  =  10|%{A'(<^)}  -  log{K{cj>)  -  A'(<^)}|  (27) 

Note  that  K{<i))  and  thus  AA(<^)  will  again  be  zero  if 

is  the  true  field  solution.  The  “dB”  error  margins  for 
the  second  and  third  order  FE/BE  method  solutions  of 
the  radiation  intensity  of  the  problem  under  considera¬ 
tion  are  shown  in  figure  9.  The  second  and  third  order 
FE/BE  method  solutions,  as  well  as  the  converged  solu¬ 
tion  of  the  radiation  intensity  are  also  shown. 

It  should  be  kept  in  mind  that  the  local  error  esti¬ 
mates,  obtained  with  equation  12,  are  dependent  on:  the 
frequency  and  polarization  of  the  incident  electromag¬ 
netic  field;  the  field  values  in  the  finite  element  under 
consideration;  the  order  of  the  approximation  functions 
in  the  finite  element  under  consideration;  the  shape  of 
the  finite  element  under  consideration;  the  accuracy  of 
the  Mrst”  FE/BE  method  solution  and  finally,  the 
accuracy  with  which  the  fields  on  the  boundary  F^^  are 
approximated  by  the  BEM  part  of  the  FE/BE  method 
solution  E^ .  The  global  error  estimate  is  dependent  on 
all  local  error  estimates  and  the  size  of  the  finite  element 
region  under  consideration. 

Bearing  all  this  in  mind  it  is  evident  that  the  local 
and  global  ERM  error  estimates  are,  in  general,  accept¬ 
ably  accurate  estimates  of  the  actual  local  and  global 
errors.  Some  of  these  above  mentioned  factors  can  be 
taken  into  account  to  improve  the  local  and  global  error 
estimates  [16,  pp.l43]  [10,  pp.l29][8] 

The  results  also  indicate  that  the  boundary  field  L^- 
norm  residual  error  estimate,  ||0pNj|^2  is  not  a  good 
quantitative  estimate  of  the  boundary  field  error.  This 
is  due  to  the  fact  that  this  error  estimate  is  a  measure  of 
the  error  in  the  boundary  field  for  the  BEM  part  of  the 
FE/BE  method  solution,  without  consideration  of  the 
actual  coupling  to  the  FEM  part  of  the  solution.  The 
actual  boundary  field  error  is,  however,  dependent  on 
the  FEM  and  BEM  part  of  the  FE/BE  method  solution. 
As  such,  the  boundary  error  estimates  must  be  used  in 
conjunction  with  the  FEM  error  estimates. 

Only  when  the  FEM  part  of  the  solution  is  of  accept¬ 


able  accuracy  and  the  BEM  part  is  highly  inaccurate  is 
this  L^-norm  error  estimate  an  accurate  quantitative  es¬ 
timates  of  the  boundary  field  error.  This  would  occur 
if  enough  finite  element  basis  functions  but  insufficient 
boundary  element  basis  functions  had  been  used  —  this 
has  not  occurred  in  the  examples  given  in  this  paper. 

It  is  evident  from  the  far-field  error  indicator  results 
(figures  6  and  9)  that,  used  together,  the  error  margins 
and  error  indicators  can  be  very  useful  in  determining 
the  accuracy  of  FE/BE  method  radar-width  or  radiation 
intensity  solutions. 


4  Conclusions 

The  variational  boundary-value  problem  formulation  of 
the  coupled  FE/BE  method  for  application  to  general 
2D  open  boundary  electromagnetic  problems  has  been 
presented.  This  formulation  has  been  developed  over 
the  past  few  decades  and  has  been  used  successfully  by 
engineers  for  numerically  solving  previously  intractable 
electromagnetic  problems.  (This  work  is  not  new  but 
was  required  as  a  basis  for  the  rest  of  the  paper.) 

ERM  local  and  global  error  estimates,  L^-norm 
boundary  field  error  estimates  and  far-field  error  indi¬ 
cators  have  been  developed  and  investigated.  This  was 
done  for  2D  FE/BE  method  solutions  of  electromagnetic 
scattering  as  well  as  radiation  problems.  The  global  EM- 
norm  ERM  error  estimate  results  obtained  for  electro¬ 
magnetic  scattering  and  radiation  problems  were  not  as 
accurate  as  the  ERM  energy  norm  error  estimates  for 
static  electric  field  problems  [19]  [16,  pp.l26].  This  is 
due  to  a  number  of  factors,  including  the  dependency 
of  the  EM-norm  error  estimates  on  the  frequency  of  the 
electromagnetic  field,  the  dependency  on  the  accuracy 
of  the  coupling  of  the  BEM  with  the  FEM,  and  the 
non-local  nature  of  electromagnetic  fields.  It  was,  how¬ 
ever,  shown  that  the  EM-norm  error  estimates  provide 
valuable  post-processed  information  regarding  the  accu¬ 
racy  and  convergence  of  the  FE/BE  method  solutions. 
This  was  shown  to  be  true  for  electromagnetic  radiat¬ 
ing  as  well  as  scattering  problems  of  arbitrary  shaped 
(lossy  and  lossless)  objects  [16].  The  local  EM-norm 
ERM  error  estimate  results  obtained  show  that  the  esti¬ 
mated  error  distributions  in  the  FEM  region  compared 
satisfactorily  with  the  actual  error  distributions.  Fur¬ 
ther  investigation  could  lead  to  improvement  of  the  lo¬ 
cal  as  well  as  global  EM-norm  ERM  error  estimates. 
(Improved  element  residual  error  estimate  methods,  for 
non-electromagnetic  problems,  have  already  been  devel¬ 
oped  [9,  20]  and  seem  to  work  well). 

Results  presented  have  show  that  the  boundary  field 
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L^-norm  error  estimate  is  not  a  quantitatively  accurate 
boundary  field  error  estimate  for  the  FE/BE  method  so¬ 
lution,  in  particular  when  the  FEM  part  of  the  solution 
is  considerably  in  error.  This  is  because  the  BEM  part 
of  the  solution  may  be  accurately  matching  an  inaccu¬ 
rate  FEM  solution  at  the  boundary.  The  boundary  error 
estimate  must  not  be  considered  in  isolation,  but  with 
due  cognisance  of  the  FEM  error  as  well. 

The  radar-width  and  radiation  intensity  error  indica¬ 
tors,  developed  for  the  FE/BE  method  solutions,  were 
used  to  obtain  “dB”  error  margins  which  proved  to  pro¬ 
vide  exceptionally  useful  post-processed  radar- width  and 
radiation  intensity  error  information.  (These  do  incor¬ 
porate,  via  the  boundary  field  derivative  error  estimates, 
information  about  the  accuracy  of  the  FEM  part  of  the 
solution  as  well) .  We  should  comment  that  Lee  in  partic¬ 
ular  has  recently  emphasized  the  role  of  “global”  errors, 
caused  by  dispersion  error:  we  have  not  considered  this 
here  [23]. 

It  was  also  shown  that  all  these  error  estimates  and 
indicators  can  be  obtained  within  negligible  computa¬ 
tional  times  compared  to  the  computational  times  of  the 
FE/BE  method  solutions.  This  is  due  to  the  nature  of 
the  error  estimate  and  error  indicator  methods  as  well 
as  the  highly  efficient  algorithms  employed  [16]. 

Adaptive  finite  element  methods  are  closely  linked  to 
a  posteriori  error  estimates  and  could  be  used  to  im¬ 
prove  the  efficiency  of  general  FEM  solutions.  The  a 
posteriori  error  estimate  methods  can  be  used  to  iden¬ 
tify  the  regions  where  the  fields  need  to  be  approxi¬ 
mated  more  accurately  (for  example  where  the  fields  vary 
more  rapidly),  and  the  finite  element  mesh  can  thus  be 
adapted  to  ensure  superior  basis  function  distributions 
in  these  regions  [10].  Although  not  shown  in  this  paper, 
the  authors  have  also  worked  on  this  topic  [16,  pp.l68]. 

In  conclusion:  the  a  posteriori  error  estimates  and  in¬ 
dicators  show  great  promise,  but  can  still  be  improved 
upon.  It  can  also  be  extended  to  a  posteriori  error  esti¬ 
mates  for  3D  FE/BE  method  solutions.  The  use  of  edge- 
based  elements  for  3D  formulations  would  not  pose  spe¬ 
cial  problems  as  far  as  error  estimation  is  concerned,  al¬ 
though  hierarchical  edge-based  elements  may  prove  more 
formidable  than  the  nodal  based  equivalents  if  the  er¬ 
ror  estimator  is  used  to  drive  an  adaptive  meshing  algo¬ 
rithm. 

Ultimately  one  would  like  to  develop  a  FE/BE  method 
solver  with  error  estimates  and  indicators  which  clearly 
identify  the  different  inaccuracies  in  the  solution  at  hand 
and  use  this  information  to  automatically  adapt  itself 
accordingly.  Such  a  solver  should  provide  accuracy  and 


convergence  information  to  the  user  providing  him  with 
a  very  useful  and  reliable  tool  for  solving  practical  elec¬ 
tromagnetic  engineering  problems.  That  this  problem 
remains  a  very  pressing  one,  and  in  particular  that  ad- 
hoc  estimates  leave  much  to  be  desired,  is  clear  from 
recent  work  [24] .  It  is  hoped  that  the  more  rigorous  esti¬ 
mate  formulations  presented  in  this  paper  will  contribute 
to  progress  in  this  field. 
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Table  1:  Percentage  error  estimates  compared  to  actual 
errors.  This  is  for  scattering  from  the  lossy  dielectric 
right-circular  cylinder  of  figure  2  excited  by  a  TEz  polar¬ 
ized,  plane  incident  field  (frequency:  2  GHz).  M  is  the 
number  of  finite  elements  and  Mb  the  number  of  bound¬ 
ary  element  used. 
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first  order  approximation  functions  and  second  or¬ 
der  approximation  functions 
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E^  second  order  approximation  functions  and  E^  third  or¬ 
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02:27:31 

00:00:40 

Table  2:  Computational  times  (on  a  HP-720  work  station)  for  the  FE/BE  method  solution  and  error  estimates  (all 
error  estimates  and  error  indicators  combined)  for  the  scattering  problem  concerning  the  right  circular  cylinder 
in  figure  2.  M  is  the  number  of  finite  element  used. 
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E^  second  order  approximation  functions  and 
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0.35 
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6.25 

5,92 
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Table  Percentage  error  estimates  compared  to  actual  errors.  This  is  for  the  electromagnetic  radiation  problem 
concerning  the  horn  antenna  of  figure  3.  Frequency:  750  MHz,  Polarization:  TEi  mode  aperture  source  field.  M 
is  the  number  of  finite  elements  and  the  number  of  boundary  element  used. 
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first  order  approximation  functions  and  second  order 
approximation  functions 


E^  second  order  approximation  functions  and  E^  third  order 
approximation  functions 


Figure  4:  Comparison  of  estimated  and  actual  local  EM-norm  error  values  for  the  FE/BE  method  solution  of  the 
electromagnetic  scattering  problem  concerning  the  right  circular  cylinder  of  figure  2.  The  incident  electromagnetic 
field  is  TE2  polarized  at  frequency  2  GHz.  M  =  501  is  the  number  of  finite  elements  used. 


Estimated  Error  Distribution 


Actual  Error  Distribution 


Figure  5:  Comparison  of  estimated  and  actual  local  EM-norm  error  distributions  for  the  FE/BE  method  solution 
of  the  fields  in  and  around  the  right  circular  cylinder  of  figure  2.  The  incident  electromagnetic  field  is  TEz 
polarized  at  frequency  2  GHz.  Dark  elements  indicate  high  error  values  and  light  elements  indicate  low  error 
values.  This  is  for  second  order  approximation  functions  and  E^  third  order  approximation  functions. 
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:  jirsi  order  approximation  functions 


second  order  approximation  functions 


Figure  6;  Comparison  of  the  radar- width  values,  obtained  analytically  and  numerically  (FE/BE  method  solution), 
for  the  scattering  problem  concerning  the  right  circular  cylinder  in  figure  2.  The  “dB”  error  margins,  A7(0),  are 
also  shown  for  first  and  second  order  FE/BE  solutions.  Frequency:  2  GHz.  Polarization:  TE^.  Number  of  finite 
elements  used:  M  =  1502.  Number  of  boundary  elements  used:  Mf,  =  90. 


first  order  approximation  functions  and  second  order 
approximation  functions 


E^  second  order  approximation  functions  and  E^  third  order 
approximation  functions 


Figure  7)  Comparison  of  estimated  and  actual  local  EM-norm  error  values  for  the  FE/BE  method  solution  of  the 
electromagnetic  radiation  problem  concerning  the  horn  antenna  of  figure  3.  Frequency:  750  MHz,  Polarization: 
TEi  mode  aperture  source  field.  M  =  450  is  the  number  of  finite  elements  used. 
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Estimated  Error  Distribution 


Actual  Error  Distribution 


Figure  8:  Comparison  of  estimated  and  actual  local  EM-norm  error  distributions  for  the  FE/BE  method  solution  of 
the  electromagnetic  radiation  problem  concerning  the  horn  antenna  of  figure  3.  Frequency:  750  MHz,  Polarization: 
TEi  mode  aperture  source  field.  Dark  elements  indicate  high  error  values  and  light  elements  indicate  low  error 
values.  This  is  for  second  order  approximation  functions  and  third  order  approximation  functions. 


Figure  9:  Comparison  of  the  radiation  intensity  values,  obtained  numerically  (FE/BE  method  solution),  for  the 
radiation  problem  concerning  the  horn  antenna  of  figure  3.  Second  and  third  order  basis  function  solutions,  the 
converged  solution  and  the  “dB”  error  margins,  are  shown.  Frequency:  750  MHz,  Polarization:  TEi  mode 

aperture  source  field.  M  and  A/5  are  the  number  of  finite  elements  and  boundary  elements  used  respectively. 
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ABSTRACT 

The  implementation  of  a  generic  mobile  telephone 
handset  model  in  the  Finite  Difference  Time  Domain 
(FDTD)  method  for  computation  of  electromagnetic  field 
distributions  is  described.  The  handset  can  be  rotated 
about  the  principal  FDTD  axes  to  achieve  any 
orientation.  The  ‘thin  wire  ’  technique  is  used  to  model 
the  antenna  in  a  way  that  ensures  that  its  correct 
electrical  length  is  maintained  despite  the  ‘staircasing’ 
effect  of  the  FDTD  grid.  Computed  predictions  jrom  the 
model  agree  with  near-field  measurements  on  a  physical 
realisation  of  the  handset,  and  accurate  near  and  far 
fields  are  calculated  for  arbitrary  orientations  of  the 
handset.  It  is  concluded  that  this  wire  rotation  technique 
has  broad  application  to  problems  involving  components 
that  are  required  to  move  in  relation  to  fixed  structures. 

1.  INTRODUCTION 

Computer  modelling  of  the  deposition  of  electromagnetic 
energy  in  the  human  head  offers  the  opportunity  to 
predict  compliance  of  radiofrequency  equipment  with 
safety  guidelines.  Datasets  representing  heads,  classified 
into  regions  with  significantly  different  dielectric 
properties,  have  been  modelled  under  irradiation  by 
plane  waves  and  dipoles  [1,2]  or  mobile  telephone 
handsets  in  simple  orientations  [3],  and  solid  spheres 
have  been  imdiated  1^  handsets  [4]  and  plane  waves 
[5,6].  Work  specifically  on  the  incorporation  of  handsets 
into  FDTD  software  has  been  performed  Luebbers  et 
al  [7],  Toftgard  et  al  [4]  and  Jensen  and  Rahmat-Samii 
[8].  With  dipole  or  handset  irradiation,  the  antenna  has 
almost  invariably  been  oriented  parallel  to  one  of  the 
principal  FDTD  axes,  so  that  its  wire  (assumed  straight) 
can  correctly  represented  in  the  FDTD  mesh. 
However,  the  standard  FDTD  methods  for  representation 
of  wires  lead  to  errors  if  the  antenna  is  not  parallel  to  a 
principal  axis  as  the  resulting  electrical  length  of  the 
?)ntenna  is  the  sum  of  a  sequence  of  orthogonal  steps 
along  FDTD  cells  (staircase  approximation).  In  contrast, 
the  physical  length  of  the  antenna  is  the  linear  distance 
between  the  initial  and  final  cells:  in  general  this  will  be 
significantly  shorter  than  the  FDTD  discretised  version. 


A  major  consequence  of  this  will  be  that  the  resonant 
frequency  of  the  antenna  will  tend  to  be  artificially 
depressed,  and  orientation-dependent,  in  the  FDTD 
representation.  The  technique  described  below 
overcomes  this  limitation,  allowing  a  handset  to  be 
located  in  a  realistic  orientation  near  a  human  head.  This 
capalnlity  is  vital  for  determination  of  the  interaction  of 
handsets  with  human  tissue  in  realistic  usage.  In 
particular,  realistic  handset  orientation  is  necessary  to 
assess  the  Specific  Absorption  Rate  (SAR)  in  the  hiunan 
head  due  to  mobile  telephones.  (It  may  be  noted  that 
Jensen  and  Rahmat-Samii  [8]  allowed  rotation  of  the 
handset  relative  to  the  head,  but  this  was  done 
rotating  the  head  rather  than  the  handset:  compensation 
for  variation  of  the  apparent  length  of  resonant 
conductors  was  thus  avoided). 

A  discussion  of  the  accuracy  that  is  estimated  for  the 
results  presented  below  is  given  in  Aj^ndix  1 . 

2.  THE  THIN  WIRE  METHOD  IN  FDTD 

The  established  method  for  representation  of  thin  wires 
in  an  FDTD  simulation  is  that  suggested  by  Holland  and 
Simpson  [9].  In  this  technique  each  element  of  the 
antenna  runs  along  an  electric  field  component  in  a  Yee 
cell  [10]  to  form  a  continuous  path  from  start  to  finish  of 
the  antenna.  At  each  time  step  the  driven  point  of  the 
antenna  is  injected  with  a  cunent,  and  the  current  on  the 
tip  of  the  antenna  is  maintained  at  zero. 

The  time-advance  equation  for  an  arbitrary  element  ‘m’ 
of  the  thin  wire,  located  in  cell  (x,y,z)=(i,j,k)  and 
pointing  along  the  x-axis,  is  given  by  [9]: 

I(m)"  =  Km)”’  + 

E,,(i,  j,  k)"-^^^  Q(m+  -Q(m)"~’^’^  ] 

where:  I(m)"  =  current  in  wire  element  m 

and  Q(m)''  =  charge  per  unit  length  on  wire  element  m 
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Ex(ij,k)“  =  X  component  of  electric  field  in  cell  (i,j,k) 
(each  of  the  above  are  for  time  instant  n) 

At  =  FDTD  time  step 

Axtw  =  ‘corrected’  length  of  thin-wire  segment  in  x- 
direction 

L  =  ‘In-cell  ’  inductance  per  unit  length  of  the  wire 

The  ‘in-ceir  inductance  per  unit  length  of  the  wire,  L, 
ensures  accurate  treatment  of  the  inductive  fields  close  to 
the  wire.  It  is  given  by  [9]: 


where  a  is  the  wire  radius,  R  is  the  effective  radius  of  an 
FDTD  cell  and  a«R. 

The  ‘corrected’  segment  length  is  found  by  multiplying 
the  standard  FDTD  cell  size  in  the  x-direction  ^  a 
correction  factor  to  reduce  the  ‘staircased’  conductor 
length  to  equal  the  correct  wire  length.  Elements  parallel 
to  the  y  or  z-axes  have  analogous  advance  equations.  The 
wire  charges  are  similarly  advanced,  using: 

(Km)""’''’  =  OCm)"-'''’  -  ( I(m)”  -  I(m-1)“)  (3) 

Axtw 

The  electric  fields  are  then  updated  by  the  currents,  after 
their  normal  FDTD  advance,  by  using: 

E,(ij,kr'''’  =  E,(ij,kr"’ .  I(mf  (4) 

e^AyAz 

where  Ay  and  Az  are  the  standard  (not  ‘corrected’) 
FDTD  cell  sizes  in  the  y  and  z  directions,  Equation  (4) 
being  for  the  case  where  the  wire  segment  is  in  the  x- 
direction. 

3.  IMPLEMENTATION  OF  ROTATABLE 
HANDSET 

To  represent  a  generic  mobile  telephone  handset,  an 
antenna  driven  as  a  monopole  and  located  on  a 
dielectric-coated  rectangular  box  was  used.  The  cell  at 
the  base  of  the  antenna  was  excited  a  current  of  the 
form  Iosin((at)  and  the  charge  in  this  cell  was  zeroed  at 
each  time  step.  The  condition  for  the  cell  at  the  tip  of  the 
antenna  was  that  the  current  be  zeroed  at  each  time  step. 

For  a  general-purpose  model  of  the  interaction  between  a 
handset  and  the  hiunan  body  it  is  normally  most 
convenient  to  keep  the  bo<^  components  fixed  in  the  co¬ 


ordinate  system  and  vary  the  orientation  of  the  handset 
to  simulate  realistic  usage.  It  is  thus  necessary  to  be  able 
to  generate  a  new,  arbitrarily  rotated,  dataset  from  an 
existing  handset  geometry  dataset.  To  define  the  handset 
structure  it  is  usually  most  convenient  to  commence  with 
the  anteima  oriented  parallel  to  an  axis  of  the  FDTD 
mesh  system  (chosen  to  be  the  x-axis  in  the  present 
work).  It  is  convenient  to  define  the  origin  to  be  at  the 
base  of  the  anterma:  orientation  is  then  defined  as  a  pair 
of  angles,  first  that  for  rotation  about  the  z-axis,  a,  and 
then  that  for  rotation  about  the  y-axis,  p  (see  Figure  1). 
This  can  readily  be  extended  to  include  an  additional 
rotation  about  the  x-axis  if  required.  The  handset 
problem  can  then  be  separated  into  two  distinct  areas: 
rotating  the  antenna  and  rotating  the  handset  body. 

3.1  Rotatable  Antenna 

To  overcome  the  problem  trf  excessive  electrical  antenna 
length  in  a  ‘staircased’  FDTD  representation,  a 
correction  factor  was  introduced,  in  which  the  values  of 
Axtw  (and  similarly  Ayt*.  and  Aztw)  used  in  Equations  (1) 
and  (3)  were  modified  from  the  basic  FDTD  cell  size  by 
multiplication  by  a  constant  fractional  factor  f*: 

Axtw  =  fw.Ax 

Aytw  =  fw.Ay 
Aztw  =  fwAz 

,  ,  Physical  distance  fromantennabase  to  tip 

where  fw=  — - - 

Staircased  distance  from  base  to  tip 

In  this  way  the  length  of  the  antenna,  as  used  in  these 
equations,  was  independent  of  rotation  angle. 

The  ‘thin  wire’  technique  requires  the  use  of  the  above 
set  of  time-advance  equations  (in  addition  to  the 
standard  FDTD  time-marching  equations)  for  each  cell 
through  which  the  wire  passes.  The  following  algorithm 
generates  the  necessary  parameters  for  use  in  Eqns  (1), 
(3)  and  (4). 

If  the  base  of  the  antenna  is  defined  as  the  origin  (0,0,0), 
the  end  point,  (Xend,yend,Zend)  is  calculated  thus: 

Xe„d=  R  cos(a)cos(p) 
ye„d=  Rsin(a) 

Zend  =  -R  cos(a)sin(P) 

where  a  is  the  angle  rotated  about  the  z-axis,  p  is  the 
angle  then  rotated  about  the  y-axis  and  R  is  the  antenna 
length.  The  FDTD  cell  indices  of  the  base  and  end  of  the 
antenna  may  then  be  ascertained. 
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It  is  next  necessary  to  calculate  the  shortest  continuous 
‘staircased’  path  of  electric  fields  between  the  two  ends, 
this  path  following  the  electric  field  components  in  the 
Cartesian  grid  of  FDTD  cells.  This  may  be  performed  by 
noting  the  indices  of  each  new  cell  traversed  by  a 
straight  line  from  base  to  tip  of  the  antenna.  The 
following  examples  show  how  the  components  and 
positions  of  electric  fields  required  to  form  such  a  thin 
wire  may  be  deduced. 

If  two  consecutive  cells  are  denoted  by  (i,j,k)  and 
(i+l,j,k),  then  the  Ex  component  located  in  cell  (i+l,j,k) 
is  required  to  link  the  two  cells.  Similarly,  the  Ex 
component  located  in  cell  (i,j,k)  is  required  to  link  the 
consecutive  cells  (i-l,j,k)  and  (ij,k)  (see  Figure  2).  In  the 
cases  where  two  indices  change  between  consecutive 
cells,  the  corresponding  two  electric  field  components 
will  be  needed.  For  example,  a  step  from  cell  (i,j,k)  to 
cell  (i+l,j+l,k)  requires  the  Ex  field  component  of  cell 
(i+l,j,k)  and  then  the  Ey  field  component  of  cell 
(i+1  j+l,k)  (see  Figure  3).  A  list  of  lattice  positions  and 
components  of  the  electric  fields  necessary  to  step  from 
base  to  tip  of  the  anterma  may  thus  be  deduced.  These 
can  then  be  used  directly  in  Equations  (1),  (3)  and  (4) 
above. 

It  was  found  that  the  formulation  of  Equations  (1)  and 
(4)  fails  if  the  chosen  direction  of  conventional  positive 
current  flow  in  the  wire  segment  under  consideration  is 
antiparallel  to  the  corresponding  unit  vector  of  the  co¬ 
ordinate  system.  To  cover  such  cases  these  equations  had 
to  be  modified  as  follows: 

I(m)"=  I(mr'  + 

x*jEx(i,j,k)"-''2  Q(m  +  l)"-*'^-Q(m)"-i^2l 

- - - ___  J 

(5) 

Ex(i,j,kr =  Ex(i,j,kr’"'  -  8(m)  I(m)"  (6) 

e^AyAz 

where  5(m)  =  +1  if  the  defined  direction  of  positive 
current  flow  in  element  ‘m’  is  parallel  to  the  unit  x- 
vector,  and  8(m)  =  -1  if  the  positive  current  flow 
direction  is  defined  to  be  antiparallel  to  the  unit  vector. 

3.2  Rotatable  Handset  Body 

As  shown  by  Luebbers  et  al  [7],  the  electrical  properties 
of  many  practical  handsets  may  be  realistically 
represented  an  antenna  on  a  metal  chassis  surrounded 
by  a  dielectric  (plastics)  case.  Consider  first  the  metal 


chassis  oriented  with  its  principal  axis  parallel  to  the  x- 
axis  as  in  Figure  1:  let  the  point  at  which  the  antenna 
and  chassis  meet  be  defined  as  position  (0,0,0)  and  let  an 
arbitrary  point  on  the  metal  chassis  be  (x„,yn„Zm)  with 
reference  to  this  origin.  As  was  the  case  for  the  antenna, 
this  point  can  be  rotated  by  a  about  the  z-axis  and  then  p 
about  the  y-axis  to  get  the  position  after  rotation 
(Xrot,yrot.Zrot)  TWs  position  is  given  by  calculating  the 
intermediate  position  (xi,yi,zi)  after  rotation  by  a  about 
the  z-axis,  then  calculating  the  position  after  a  further 
rotation  by  p  about  the  y-axis,  viz. : 

xi  =  ro  cos(0o  +  a) 
yi  =  ro  sin(0o  +  a) 

where  0o  =  tan  ' (yn,/Xn,) 
and  ro  =  (x„^-i-ymY’ 

Then: 

x„t  =  risin(0i  +  p) 
yrot  =  yi 

z„t  =  ricos(0i  +  P) 

where  0i  =  tan''(zi/xi) 
and  ri  =  (xi^  +  zi  Y 

Each  point  represents  a  location  on  the  metal  chassis 
after  rotation.  The  lattice  cell  indicated  Ity  each  such 
point  is  easily  calculated  and  its  position  must  be  stored 
for  later  use  during  time  stepping.  At  each  time  step,  the 
tangential  electric  field  components  on  the  metal  surface 
must  be  zeroed.  In  a  staircased  representation  this  is 
achieved  by  zeroing  all  three  orthogonal  components  in 
any  cell  deemed  to  contain  metal. 

Similarly,  the  dielectric  case  of  the  handset  may  be 
rotated.  In  this  case  it  is  necessary  only  to  set  correct 
electrical  properties  for  the  indicated  lattice  cells.  In 
general,  any  dielectric  or  perfectly  conducting  metal 
shape  may  be  rotated  in  this  manner. 

4.  RESULTS 

Numerical  tests  were  undertaken  with  an  FDTD  program 
based  on  THREDE  [11],  but  using  Fang-Mei 
superabsorbing  bormdary  conditions  as  the  interface  to 
free  space  [12].  In  addition,  a  near-to-far  field 
transformation  algorithm  (see  below)  and  the 
modifications  described  above  were  incorporated  into  the 
program.  The  computer  was  a  Sun  Sparcstation  10  with 
96Mbyte  main  memory. 
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4.1  Near  Fields 

A  batteiy-powered  transmitter  simulating  the  generic 
shape  of  typical  commercial  handsets  was  constructed. 
This  radiated  1.8GHz  CW  from  a  ^4  monopole  located 
on  the  top  surface  of  a  rectangular  case.  The  case  was 
made  of  plastics  material  (with  Er=2.0,  as  in  [7])  with  a 
conducting  inner  coating.  This  was  placed  in  an 
anechoic  chamber,  and  near-field  measurements  of  the  x- 
component  of  the  electric  field  were  made  over  the 
volume  shown  in  Figure  4.  An  FDTD  representation  of  a 
handset  with  the  same  dimensions  was  constructed  and 
the  same  field  component  was  calculated.  Figure  5  shows 
the  electric  field  intensity  in  a  series  of  planes  at  regular 
intervals  away  from  the  handset.  The  measured  results 
can  be  seen  to  compare  very  closely  with  results  from 
FDTD  analysis  in  the  region  close  to  the  handset.  The 
discrepancy  appears  to  increase  somewhat  at  larger 
distances  (maximum  error  within  3dB),  but  the 
agreement  can  still  be  said  to  be  good,  given  the 
difficulties  of  representing  the  exact  structure  of  a 
simulated  handset  assembled  from  commercial 
components,  and  errors  inherent  in  field  strength 
measurements  with  a  small  dipole  probe.  The  results  can 
be  said  to  demonstrate  the  validity  of  Luebbers’  model  of 
an  antenna  on  a  plastic  coated  metal  chassis  [7]  in  the 
near  field. 

Figure  6  shows  the  total  electric  field  magnitude,  as 
calculated  1^  the  FDTD  method,  in  the  arrangement 
described  above.  Also  shown  are  the  results  calculated  by 
rotating  the  FDTD  handset  model  by  45°  about  the  z- 
axis.  The  results  for  the  rotated  handset  are  seen  to  give 
excellent  agreement  with  those  for  the  axially-aligned 
handset. 

The  far  field  was  also  calculated  for  the  axially  aligned 
and  rotated  handsets  using  the  near-to-far  field 
transformation  method  described  by  Taflove  & 
Umashankar  [13].  Figure  7  shows  the  far-field  radiation 
pattern  in  the  x-y  plane  for  both  orientations:  the  results 
from  the  rotated  handset  have  been  rotated  back  by  45° 
to  enable  direct  comparison  to  be  made.  The  amplitudes 
are  in  good  agreement,  with  a  slight  angular  offset. 
Empirically,  the  angle  of  error  has  been  found  to  be  of 
the  order  of  1/N  radians,  where  N  is  the  munber  of  cells 
in  the  FDTD  antenna.  This  indicates  that  quantization  is 
the  cause  of  the  angular  error;  use  of  finer  FDTD  meshes 
reduces  it. 

4.2  Output  Power 

The  power  radiated  the  FDTD  handset  model  was 
measured  for  a  range  of  orientation  angles.  Although 


the  radiation  patterns  were  in  excellent  agreement  for  all 
orientation  angles,  it  was  foxmd  that  the  actual  power 
radiated  was  dependent  upon  orientation.  Let  the  angle  p 
now  be  defined  as  the  angle  between  the  antenna  and  the 
x-axis  in  the  x-z  plane  (see  Figure  1).  The  FDTD 
handset  model  was  set  to  a  number  of  values  of  p 
between  0  and  Till  and  the  radiated  power  was  computed 
in  each  case.  The  power  calculation  was  performed  by 
summing  the  outgoing  component  of  the  Poynting  vector 
over  a  closed  siuface  surrounding  the  handset.  Figure  8 
shows  the  power  transmitted  as  a  fimction  of  angle  when 
the  monopole  was  fed  with  lA  peak  CW  at  1.8GHz.  The 
FDTD  cells  were  cubes  with  side  lengths  of  2.5mm, 
corresponding  to  approximately  >733  in  free  space. 

The  radiated  power  can  be  seen  to  vary  by  around  20% 
as  the  angle  is  varied.  The  power  output  is  symmetrical 
about  the  7t/4  orientation  as  would  be  expected  in  a  cubic 
FDTD  mesh.  It  was  suspected  that  the  variation  might  be 
due  to  quantization  errors  in  the  length  of  the  antenna 
and  the  handset  Ixxfy.  The  anteima  is  required  to  be 
terminated  by  zeroing  the  current  in  its  end  segment.  As 
different  orientations  imply  differing  values  of  Axtw,  Ay^^ 
and  Aztw  the  ‘thin  wire’  antenna  may  be  expected  to  be  a 
source  of  error,  since  end  segments  of  varying  lengths 
would  be  partially  excluded  from  the  electrical  length. 

Fig.  9  shows  the  currents  along  the  antenna  for 
orientations  between  0  and  45°.  The  results  are 
symmetrical  about  45°  so  results  from  45  to  90°  are  not 
shown.  The  currents  follow  a  similar  pattern,  but  with 
discrepancies  of  up  to  10%,  which  might  account  for  the 
20%  variation  in  total  power.  To  stucfy  this  phenomenon 
fiirther,  the  same  problem  was  re-computed  using  a 
reduced  cell  size  of  1.25mm  (A,/67).  Figures  10  and  11 
show  the  power  and  current  for  a  range  of  orientation 
angles  at  this  resolution.  The  variation  in  power  (Figure 
10)  is  reduced  to  6%  maximiun,  and  the  variation  in 
current  (Figure  1 1)  is  nowhere  higher  than  7%  along  the 
lengths  of  the  wires.  The  relatively  large  drop  in  power 
variation,  compared  with  the  result  for  Ay33  resolution, 
suggests  that  the  handset  Ixxfy  may  be  equally 
responsible  for  power  variations.  The  finer  cell  size 
enables  the  representation  of  a  rotated  body  to  be 
electrically  more  exact.  The  marked  improvement  in 
correlation  of  both  power  outputs  and  currents,  as 
resolution  improves,  indicates  that  the  small 
discrepancies  originate  from  quantization  errors,  not 
errors  in  the  method.  It  is  important  to  note  that,  when 
field  patterns  are  normalised  to  1  Watt  output  (as  is 
normal  in  the  Specific  Absorption  Rate  (SAR) 
calculations,  for  which  the  rotatable  handset  model  has 
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been  developed),  both  near  and  far-field  patterns  are  in 
excellent  agreement  at  all  orientation  angles. 

5.  CONCLUSION 

It  has  been  shown  that  the  ‘thin  wire’  technique  can  be 
used  as  a  basis  for  modelling  an  arbitrarily-oriented 
handset  (or  some  other  structure  made  of  dielectrics  and 
metal)  in  an  FDTD  lattice.  Both  near  and  far  fields  are  in 
good  agreement  with  the  results  for  a  handset  in  the 
‘conventional’  orientation,  i.e.  with  principal  dimensions 
parallel  to  the  Cartesian  axes,  and  the  near  fields  also 
agreed  well  with  physical  measurements.  This  technique 
has  particular  application  to  the  calculation  of  ^e 
interaction  between  realistically-oriented  portable 
handsets  and  the  hiunan  bo<^.  The  technique  for  rotating 
a  solid  object  in  an  FDTD  lattice  has  broad  application  to 
problems  involving  two  or  more  physical  entities  that 
may  be  re-positioned  in  relation  to  each  other  (e.g. 
helicopter  rotors  and  microwave  oven  interiors).  One  or 
more  of  the  objects  can  be  positioned,  and  repositioned, 
relative  to  others  without  significant  restructuring  of  the 
program  or  input  data. 
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APPENDIX  1 

Quantitative  Accuracy  Statements 

1.  An  estimate  of  the  accuracy  of  the  results  presented 
here  can  be  made  from  examination  of  the  comparisons 
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presented  in  Figs  5,  6,  7,  9  and  11,  and  from  a 
comparison  of  Figs  8  and  10. 

In  Fig.  5,  it  is  seen  that  the  shapes  of  the  computed  near- 
field  distributions  agree  well  with  the  measured  results 
(normalised  to  equality  at  identical  central  points)  and 
the  maximum  magnitude  discrepancy  is  within  3dB. 
E?^rience  suggests  that  this  discrepancy  is  partly 
caused  by  the  limitations  of  the  quantised  FDTD 
representation  of  the  handset  structure  and  partly  by 
measiu'ement  uncertainty.  Fig.  6  indicates  that  computed 
near  field  results  are  stable  to  within  IdB  imder 
rotational  transformation  with  respect  to  the  FDTD 
mesh. 

Fig.  7  indicates  stability  in  the  peak  far  field  gain  of 
better  than  0.5dB  under  rotational  transformation.  Errors 
of  the  order  of  lOdB  may  occur  in  deep  nulls,  but  this  is 
a  common  problem  which  could  be  reduced  using  a 
smaller  cell  size. 

Comparison  of  Figs  8  and  10  shows  that  the  cell  size  of 
Ay33  used  to  produce  the  above  results  also  leads  to 
variations  of  about  ±0.5dB  in  the  radiated  power  under 
rotational  transformation,  but  that  halving  the  cell  size 
reduces  this  to  ±0. 12dB. 

Fig.  9  indicates  maximum  variations  of  about  ±0.2dB  in 
the  antenna  current  distribution  under  rotational 
transformation  with  a  cell  size  of  >./33,  reducing  to 
±0.13dB  with  X/67  cells  (Fig.  11). 


Fig.l  The  simulated  generic  handset,  associated  co¬ 
ordinate  system  and  rotation  angles. 


2.  These  estimates  are  based  on  measurement  (Fig.  5), 
co-ordinate  transformation  and  variation  of  discretisation 
intervals. 

3.  The  kinds  of  problems  for  which  the  above 
error/accuracy  statements  can  be  made  for  this  program 
include  the  modelling  of  cunent  distributions,  near  and 
far  fields  for  arbitrary  structures  of  conductor  and 
dielectric,  with  maximum  dimensions  up  to  about  a 
wavelength. 

4.  Nominal  sampling  densities  required  to  achieve  these 
estimated  accuracies  are  Xy33  in  most  cases  (X/67  where 
indicated). 

5.  The  operation  count  needed  to  exercise  the  model 
reported  here  is  estimated  nominally  to  be  of  the  order  of 
10'¥*,  where  f  is  in  Hz. 

6.  The  variable  storage  needed  to  exercise  this  model  is 
estimated  nominally  to  be  IMword  for  X/33  sampling 
density  (8Mwords  for  X/67).  This  is  based  on  the 
modelling  of  a  simulated  handset  at  1.8GHz,  where  ^33 
=  5mm.  The  model,  with  the  same  absolute  dimensions 
and  sampling  interval,  could  be  used  down  to  about 
IGHz,  below  which  the  free-space  margin  between 
conductors  and  the  absorbing  boundary  would  have  to  be 
increased;  it  should  also  remain  stable  up  to  about  4GHz 
(based  on  sampling  interval  considerations),  although 
substantial  errors  would  begin  to  appear  in  some 
parameters. 


1 

Ex 

s/ 

Ex(i+l^jyk) 

cell 

? 

cell  i,j,k 

cell  i+lj,k 

Fig.2  Rectilinear  linkage  of  electric  field  components 
between  adjacent  FDTD  cells:  a  continuous 
sequence  of  such  links  is  needed  to  r^resent  a 
thin  wire. 
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Ey  (i+Ej+Ek) 


Fig.  3  Angled  linkage  of  electric  field  components  between  non-rectilinear  sequences  of  FDTD  cells. 


Fig.4  Dimensions  of  handset  and  volume  over  which  fields  were  computed  and  measured.  Dimensions  are  given  in 
ffee-space  wavelengths  at  1.8GHz. 
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Time  domain 


Measured 


Fig.  5  Magnitude  in  dB  of  the  x-component  of  the  electric  field  over  several  planes  normal  to  the  z-axis:  (a)  as 

computed  by  FDTD;  (b)  as  measured.  Both  plots  are  normalised  to  their  respective  amplitudes  at  the  centre  of 
the  plane  0. 17m  from  the  handset. 


Distance  from  Handset 


Magnitude  (in  un-normalised  dB)  of  the  total  electric  field  (|E|)  over  several  planes  normal  to  the  z-axis  as 
computed  by  FDTD;  (a)  handset  axially-aligned;  (b)  handset  rotated  45°  about  z-axis. 


tUA  rotated 


Fig.7  Far-field  pattern  in  plane  z=0,  as  computed  by  Fig.8  Total  radiated  power  predicted  by  FDTD  (mesh 
FDTD  with  handset  axially  oriented  (solid  interval  =  X/33),  for  constant  antenna  feec^wint 

curve)  and  rotated  45°  about  z-axis  (the  broken  current,  as  a  function  of  handset  rotation  angle, 

curve  shows  the  result  as  rotated  back  ly  -45°  to 
facilitate  comparison).  Scale  is  in  dB  relative  to 
the  same  reference  level,  with  a  constant  offset 
added  to  make  all  results  positive. 


Fig,  9  Variation  of  computed  current  distribution  in  the  monopole  with  rotation  of  the  handset,  for  FDTD  mesh 
interval  of  X/33.  Normalised  current  is  displayed  as  the  ratio  of  local  r.m.s.  value  to  peak  driving  current  at 
base  of  antenna. 
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Fig.  10  Total  radiated  power  as  a  function  of  orientation  angle,  for  FDTD  mesh  interval  of  X/67. 


Fig.  1 1  Variation  of  computed  current  distribution  in  the  monopole  with  rotation  of  the  handset,  for  FDTD  mesh 
interval  of  ^767  (Current  normalised  as  in  Fig.  9). 
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Abstract 

A  Magnetic  Field  Iterative  Technique  (MFIT)  is  developed 
as  a  technique  for  improving  the  results  from  high  frequency 
based  prediction  methods.  The  technique  combines  the 
accuracy  of  low  frequency  methods  with  the  speed  of  high 
frequency  methods  to  develop  a  contraction  mapping  based, 
iterative  solver,  which  is  directly  parallelizable  on  current 
massively  parallel  processing  machines.  In  this  paper,  the 
MFIT  is  developed  from  a  Magnetic  Field  Integral  Equation 
(MFDE)  and  its  contraction  mapping  properties  are 
discussed.  Results  obtained  using  the  technique  and 
illustrations  of  the  technique’s  convergence  properties  are 
presented  for  both  two  and  three  dimensional  perfect 
electrically  conducting  (PEC)  targets. 

Introduction 

High  frequency  based  numerical  prediction  methods  such  as 
Physical  Optics  -  Physical  Theory  of  Diffraction  (PO-PTD) 
[1]  and  Geometrical  Optics  -  Geometrical  Theory  of 
Diffraction  (GO-GTD)  [2]  have  been  shown  to  be  successful 
at  approximating  the  induced  current  distributions  on 
electrically  large  targets  (e.g.,  three  dimensional  bodies  with 
at  least  one  dimension  greater  than  approximately  ten 
wavelengths).  For  these  electrically  large  bodies,  the  high 
frequency  series  expansion  associated  with  their  scattering 
properties  converges  quickly  and  scattering  centers  located 
about  the  bodies  can  be  easily  identified.  However,  as  the 
electrical  size  of  these  perfect  electrically  conducting  (PEC) 
bodies  begin  to  decrease,  the  ability  to  clearly  define  the 
required  scattering  centers  becomes  difficult  and  the 
convergence  rate  of  the  associated  series  expansions 
decreases  rapidly.  Most  often,  more  exact  numerical 
methods  such  as  the  Method  of  Moments  (MoM)  [3,4]  and 
the  Finite  Element  Method  (FEM)  [5]  have  been  utilized  for 
these  particular  situations.  While  these  methods  (MoM  and 
FEM)  have  been  proven  to  be  both  accurate  and  reliable  for 


PEC  bodies  which  are  less  than  approximately  three 
wavelengths  on  a  side,  these  methods  become  difficult  to 
implement  for  much  larger  bodies.  Large  memory 
requirements  coupled  with  long  execution  times,  make  the 
utilization  of  these  methods  difficult  for  large  bodies. 

The  Magnetic  Field  Iterative  Technique  (MFIT)  is 
developed  here  for  those  PEC  bodies  whose  electrical  size 
falls  into  the  awkward  gap  which  exists  between  the  two 
groups  of  numerical  methods.  Currently,  there  are  few 
options  that  exist  for  efficiently  calculating  the  induced 
current  distributions  on  these  particular  PEC  bodies.  The 
following  method  is  a  hybrid  numerical  technique  which 
takes  advantage  of  the  desired  features  from  both  the  high 
frequency  and  nearly  exact  solution  methods.  The  method  is 
an  extension  of  previous  research  performed  by  Thiele  et  al. 
[6-8]  on  hybrid  and  iterative  methods  for  solving  electrically 
large  problems. 

Mathematical  Development 

For  closed  volumetric  PEC  bodies,  a  Magnetic  Field 
Integral  Equation  (MFIE)  may  be  stated  in  terms  of  the  total 
magnetic  field  as  [3]: 

W{7)  =—UhxH\r'))^'G{r,r')ds+H'{r)  (1) 

where  lim —  [ (n  X  (r'))  X  V'G{r  ,r')ds 

An ' 

=  |(F^(r)-«(nF^(r))  (2) 

S  is  the  closed  volumetric  PEC  surface,  G  is  the  standard 
free  space  Green’s  function,  T  is  the  total  field,  and  I  is  the 
incident  field.  Equations  1  and  2  are  in  the  form  of  a 
contraction  mapping  [9,10],  and  it  may  be  shown  that  direct 
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incident  field.  Equations  1  and  2  are  in  the  form  of  a 
contraction  mapping  [9,10],  and  it  may  be  shown  that  direct 
iteration  of  these  equations  will  yield  monotonic  mean- 
square  convergence  provided  that  the  integral  operator  is 
bounded  by  unity.  If  the  integral  operator  does  satisfy  this 
bound,  then  the  above  convergence  is  guaranteed 
independent  of  the  initial  guess.  While  this  unity  bound 
may  be  shown  for  explicit  geometries,  a  general  proof  does 
not  yet  exist.  Extensive  discussion  regarding  the  guarantee 
of  convergence  may  be  found  in  references  [11-14]. 

Since  equations  1  and  2  are  normally  implemented  in 
discretized  form,  it  is  best  to  discuss  the  issue  of 
convergence  in  terms  of  the  discretized  equations. 
Applying  a  pulse-basis  point-matching  expansion  to  the 
above  MFIE  results  in  the  following  discretized  equation 
where  delta  is  the  area  of  the  discretized  facets  (see  Figure 
1): 


H:  =Y.^„GiR„„) 


J.,  +-J...+HL 


(3) 


where  = - irz — ; - 

i«..i 

Expressing  the  discretized  MFIE  in  matrix  form  transforms 
the  problem  of  establishing  convergence  from  bounding  the 
integral  operator  to  determining  the  spectral  radius  of  the 
reaction  matrix  [Zj^,]  (Equation  6).  If  the  spectral  radius  of 
the  reaction  matrix  is  less  than  unity,  then  the  discretized 
MFIE  is  a  true  contraction  mapping  and  direct  iteration  of 
the  discretized  equation  will  guarantee  monotonic  mean 
square  convergence. 


(4) 
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For  cases  where  the  spectral  radius  is  not  less  than  unity,  it 
is  possible  to  shift  the  spectral  radius  of  the  matrix  [Z„„J  by 
using  a  relaxation  parameter  [1 1-14].  The  implementation 
of  the  relaxation  parameter  a  is  given  in  Equation  7  and  its 
numerical  value  is  a  function  of  the  largest  eigenvalue  of  the 
[Zn^J  matrix.  Thus,  for  all  given  reaction  matrices,  it  is 
possible  to  obtain  the  desired  solution  through  direct 
iteration  with  the  guarantee  of  monotonic  mean-square 
convergence  if  a  is  properly  chosen  as  required  by  the 


following  equation  where  L  is  the  linear  operator  used  in 
Equation  6: 

)  +  (1  - G)  )//*  0  <  1  (7) 

While  the  guarantee  of  convergence  is  independent  of  the 
initial  guess,  the  rate  of  convergence  is  not;  and  thus,  for  the 
MFIT  to  serve  as  an  efficient  method  to  obtain  the  induced 
current  distributions  on  volumetric  PEC  bodies,  the  iterative 
procedure  must  be  initiated  with  the  best  guess  available 
such  that  the  computation  time  associated  with  the 
calculation  of  the  initial  guess  is  small  compared  to  the 
computation  time  associated  with  the  iterative  procedure. 
For  this  reason,  high  frequency  techniques  such  as  GO-GTD 
and  PO-PTD  are  used  to  generate  the  initial  guess.  Both  of 
these  techniques  have  negligible  computer  run  times  and  are 
capable  of  producing  approximations  for  the  induced 
current  distributions  which  are  in  the  realm  of  rapid 
convergence.  It  is  the  use  of  the  high  frequency  prediction 
methods’  results,  as  an  initial  guess  to  the  iterative  process, 
which  creates  the  hybrid  method.  By  initiating  the  iterative 
process  with  these  results,  the  large  amount  of 
computational  work  associated  with  calculating  the 
principle  current  distributions  using  the  more  exact  methods 
is  removed  from  the  solution  process,  leaving  the  iterative 
process  to  fine  tune  the  results  rather  than  regenerating  the 
first  order  solutions  which  constantly  reoccur.  In  essence,  a 
high  frequency  prediction  method  is  used  to  quickly  obtain 
the  general  portion  of  the  solution  and  the  iterative  solution 
method  is  used  to  fine  tune  the  general  solution  to  the 
particular  PEC  body. 


m*  facet 


Figure  1  Discretized  Geometry 

2-Diniensional  Example 

For  the  2-dimensional  transverse  electric  (TE)  case,  the 
previous  vector  MFIE  (Equations  1  and  2)  simplifies  to  the 
following  scalar  integral  equation  [6]. 
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Hl{f)  =  \Hl{r)G{r,r')cos{Q)dUHi (r)  (8) 

Applying  a  pulse-basis  point-matching  expansion  to  the 
above  integral  operation  results  in  the  following  series 
expansion  [3]. 

hI  =  -^Y.^Xg{\r4)cos{q,„„)+1^h:^hI  (9) 

This  2-dimensional  series  expansion  was  applied  to  a  1 
wavelength  diameter  PEC  cylinder  (see  Figure  2).  The  2- 
dimensional  contour  was  discretized  at  twenty  points  per 
wavelength,  and  the  iterative  procedure  was  initiated  with 
its  high  frequency  Physical  Optics  approximation  [1]. 
Figure  4  shows  the  converging  currents  for  the  circular 
cylinder  starting  with  the  initial  guess,  then  the  current  after 
1,  5,  and  10  iterations  ending  with  the  “converged”  current. 
For  this  example,  convergence  was  based  upon  a  tight 
mean-square  error  norm  of  0.001.  The  converged  solution 
was  obtained  after  19  iterations  and  the  vast  improvement 
after  only  5  iterations  should  be  noted.  As  can  be  seen,  the 
monotonic  mean  square  convergence  associated  with 
contraction  mappings  is  present  in  the  converging  iterative 
solution. 


Figure  2  One  Wavelength  PEC  Cylinder  Geometry 

3-Dimensional  Example 

A  1  wavelength  PEC  cube,  with  a  normally  incident 
transverse  electric  (TE)  plane  wave  is  shown  as  a  3- 


dimensional  example  (see  Figure  3).  The  cube  geometry 
was  selected  because  of  the  non-uniform  and  singular 
current  distributions  that  exist  about  its  sharp  edged 
structure.  These  particular  current  distributions  have  been 
traditionally  difficult  to  calculate  and  would  serve  as  a 
rigorous  test  on  the  contraction  mapping  properties  of  the 
MFIE  formulation.  The  cube  surface  was  facetized  with 
approximately  1,300  triangular  facets  and  Xpatch  [15],  a 
high  frequency  ray  based  prediction  code,  was  used  to 
provide  the  initial  current  distribution  (see  Figure  5a).  As 
can  be  seen  from  Figure  5a,  while  the  high  frequency 
prediction  does  a  good  job  of  predicting  the  current 
distribution  on  the  front  face  of  the  cube,  it  fails  to  predict 
any  current  distributions  on  the  remaining  five  faces  of  the 
cube. 

Performing  a  single  iteration  on  the  high  frequency 
predicted  current  distribution  (Figure  5b)  results  in  a  non¬ 
zero  current  distribution  on  all  faces  of  the  cube  surface. 
Further  iterations  (Figures  5c  and  5d)  again  illustrate  the 
contraction  like  mapping  properties  of  the  MFIE 
formulation  and  the  convergence  of  the  induced  current 
distributions  along  the  surfaces  of  the  PEC  cube. 
Comparing  the  high  frequency  predicted  current  distribution 
(Figure  5a)  with  the  converged  current  distribution  obtained 
after  20  iterations  (Figure  5d)  shows  a  substantial  difference 
in  all  but  the  front  face  of  the  cube.  Finally,  it  should  be 
noted  that  all  of  the  non-uniform  and  singular  current 
distributions  are  present  and  observable  in  the  converged 
current  distribution,  and  that  these  difficult  current 
distributions  had  no  detrimental  effect  on  the  convergence 
of  the  iterative  procedure. 


d=U 


Figure  3  One  Wavelength  PEC  Cube  Geometry 


68 


(c)  10  Iterations  (c)  20  Iterations 

Figures  5a-d  Current  Magnitudes  on  a  1  Wavelength  Cube 


7  0 


Conclusions  and  Comments 


References 


Based  upon  the  results  presented,  the  MFIT  shows  promise 
as  an  efficient  method  for  bridging  the  gap  between 
traditional  high  frequency  and  low  frequency  numerical 
techniques.  Although  this  hybrid  method  is  designed  for 
bodies  which  are  3  to  10  wavelengths  in  size,  the  1 
wavelength  bodies  selected  in  the  examples  were  chosen 
because  of  their  ability  to  compactly  illustrate  the 
fundamental  principles  underlying  the  MFIT.  The  use  of 
high  frequency  methods  to  provide  an  initial  guess  to  an 
iterative  solution  formulation  of  low  frequency  methods 
allows  the  combination  of  the  best  of  both  techniques.  The 
information  readily  obtained  using  high  frequency  methods 
greatly  reduces  the  run  time  associated  with  the  low 
frequency  methods.  Hence,  the  low  frequency  portion  of 
this  hybrid  method  is  used  solely  to  obtain  the  desired 
solution  accuracy  and  valuable  computer  time  is  not  wasted 
on  obtaining  portions  of  the  solution  which  are  well  known 
or  more  readily  obtainable.  Furthermore,  iterative  solvers 
for  systems  of  linear  equations  are  highly  parallelizable  on 
current  massively  parallel  processing  machines,  thus  further 
reducing  the  run  times  associated  with  the  low  frequency 
portion  of  the  hybrid  technique. 

The  contraction  mapping  principle,  associated  with  the 
formulation  of  the  MFIT  presented,  removes  many  of  the 
dangers  commonly  encountered  when  working  with 
iterative  solvers  [9-10].  The  use  of  the  relaxation  parameter 
a  allows  generic  reaction  matrices  to  be  quickly 
reformulated  into  contraction  mapping  based  systems.  For 
more  information  on  the  use  of  relaxation  parameters  to 
ensure  contraction  mapping  properties,  the  authors 
recommend  references  [11-14]  and  the  references  therein. 

In  summary,  a  procedure  has  been  presented  which  uses  a 
magnetic  field  integral  equation  in  an  iterative  process.  The 
MFIE  is  discretized  and  a  relaxation  parameter  is  used  to 
ensure  the  spectral  radius  of  the  reaction  matrix  is  less  than 
unity.  The  discretized  MFIE  has  contraction  mapping 
properties  such  as  guaranteed  monotonic  mean  square 
convergence.  The  iterative  process  starts  with  a  first  order 
approximation  to  the  current  by  using  a  high  frequency 
technique  such  as  physical  optics  or  geometrical  optics 
(shooting  and  bouncing  rays)  rather  than  starting  at  zero 
current.  The  use  of  a  high  frequency  technique  as  a  starting 
point  in  the  iterative  process  significantly  speeds  up  the 
iterative  convergence  but  is  not  necessary  to  obtain  a 
converged  result  since  that  is  guaranteed  by  the  contraction 
mapping.  A  future  paper  will  deal  with  two  MFIE 
formulations  in  the  MFIT  and  explore  their  convergence 
properties  in  more  detail. 
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Abstract- An  integral  formulation  is  presented  and  used  for 
the  analysis  of  the  Remote  Field  Eddy  Current  (RFEC)  tran¬ 
sient  phenomena.  The  model  presented  allows  simple  and 
accurate  modeling  of  RFEC  systems  in  the  presence  of 
axisymmetric  defects. 

The  method  has  been  validated  flrst  by  comparison  of  the 
computed  results  with  results  obtained  with  an  analytical 
model.  Then  the  method  was  used  to  show  that  the  Remote 
Field  effect  is  mainly  due  to  a  direct  induction  effect  rather 
than  propagation  effects. 

The  Remote  Field  effect  phenomena  are  usually  analysed 
using  a  steady  state  domain.  However,  the  results  obtained  here 
have  shown  that  detection  of  defects  is  not  due  to  a  steady  state 
phenomenon,  but  depends  mainly  on  the  transient  signal. 


defects  on  either  surface  of  the  tubes.  This  model  is  not 
supported  by  experimental  data  or  by  the  basic 
electromagnetic  theory  for  low  frequency  fields.  Existing 
models  for  the  far  field  effect  assume  that  all  testing  is  done 
at  such  low  speeds  that  the  motion  of  the  probe  does  not 
affect  the  signal  obtained.  This  is  far  fi-om  being  realistic. 
At  low  frequencies  the  speed  effects  on  the  signal  are  large 
and,  in  some  cases,  may  even  dominate.  In  this  paper,  we 
propose  an  integral  model  that  incorporates  speed  effects,  to 
correctly  analyse  the  influence  of  speed  on  the  calculated 
data. 

2.  MODEL 


1.  INTRODUCTION 

The  remote  field  effect  is  a  Non  Destructive  Test  method 
used  for  the  testing  of  thick  ferromagnetic  materials  for 
deep  defects  and  in  particular  corrosion  effects  on  the  outer 
surface  of  tubes.  This  is  particularly  useful  for  gas  lines  but 
also  for  water  lines  and  other  thick  walled  tubes  like  oil  well 
casings.  The  method  relies  on  the  use  of  veiy  low  intensity, 
low  frequency  electromagnetic  fields  which,  unlike  other 
types  of  tests  are  equally  sensitive  to  defects  on  the  outer 
surface  of  the  tubes  as  on  the  inner  surface.  Thus,  testing  of 
thick  materials  can  be  easily  performed  from  the  accessible 
inner  surface  [1].  The  RFEC  testing  probe  is  made  of  two 
coils  separated  at  a  distance  of  about  two  coil  diameters  and 
moving  as  a  unit. 
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Fig.  1  RFEC  probe  configuration 


One  coil  (the  exciting  coil)  is  driven  with  a  sinusoidal 
current,  and  induces  an  output  signal  in  the  second  coil  (the 
pick-up  coil).  The  probes  move  inside  the  tube  at  a  constant 
velocity  and  the  exciting  current  frequency  is  normally  be¬ 
low  100  Hz.  The  theoretical  basis  of  the  remote  field  effect 
is  not  well  established.  In  particular,  some  models  [2-4]  call 
for  propagation  effects  to  e.\plain  the  equal  sensitiviri'  to 


Due  to  the  axial  symmetry  of  the  system,  we  divide  the 
tube  into  elementary  rings  in  which  a  uniform  current  is  as¬ 
sumed. 


Fig  2  Subdivision  of  the  tube  into  concentric  rings 


Wc  write  the  Ohm's  law  at  a  point  inside  a  generic  ring 


as; 


a  dt 

and  then  integrate  both  sides  of  this  equation  over  the  vol¬ 
ume  of  the  ring  to  obtain  the  integral  equation  [5]; 

,  .  die  ^  dMpe  ,  dik  dij  ^  rry. 

dt  dx  fc=\  dt  dt 

This  can  be  viewed  as  the  electric  equilibrium  equation  of 
the  equivalent  circuital  loop  made  of  a  resistance  R,  a  self 
inductance  L  due  to  the  ring  itself  and  the  mutual  induc¬ 
tances  M  between  the  ring  and  all  other  rings  that  make  up 
the  tube  and  the  coils.  In  Eq.  (2)  the  subscripts  k  identifies 
the  rings  I  through  n,  e  identifies  the  exciting  coil,  \  the 
current  in  the  k-th  loop  of  the  network  and  the  current  in 
the  exciting  coil. 

To  calculate  the  signal  we  consider  only  the  n  significant 
rings  close  to  the  coils.  Writing  this  equation  inside  everj^ 
ring,  we  obtain  a  system  of  n  equations,  that  can  be  viewed 
as  the  equilibrium  equations  of  the  equivalent  network 
shown  in  fig.  3.  The  equations  thus  obtained  can  be  easily 
integrated  using  a  classical  Runge-Kutta  method. 
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Fig.  3  Equivalent  electric  network 

Axisymmetric  defects  in  the  tube  are  modeled  by  simply 
removing  the  rings  corresponding  to  the  volume  of  the  de¬ 
fects.  The  corresponding  loops  in  the  equivalent  network  are 
also  removed.  The  output  signal  Vp  is: 

f,  dlk  dMpk  ..  die 

k=i  dx  dt  (3) 

The  currents  1^  are  found  from  the  solution  of  the  sys¬ 
tem  of  equations  for  the  equivalent  network  in  fig.  3. 


3. 


RESULTS 


considerable  influence  on  the  correctness  of  the  calculated 
data.  To  obtain  good  results,  the  distance  between  the  excit¬ 
ing  coil  and  the  edge  of  the  tube  must  be  at  least  twice  the 
distance  between  the  exciting  and  pick-up  coils.  Next,  we 
analyze  the  effect  of  discretization  on  the  calculated  data 
using  a  400mm  tube  and  three  levels  of  discretization.  Fig.  6 
shows  that  the  calculated  results  agree  well  with  the  ana¬ 
lytical  results  even  for  a  coarse  discretization  of  50  x  3  rings. 


Fig.  6  Induced  voltage  vs.  coil  velocity,  tube  without  defect 


The  model  presented  above  has  been  tested  by  analyzing 
the  geometry  of  TEAM  Problem  9  [6],  shown  in  fig.  4.  The 
results  for  the  tube  without  defects  were  compared  with  re¬ 
sults  obtained  from  an  analytical  solution. 
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Fig.  4  Geometry  of  the  Team  problem  9 

We  first  analyze  the  behaviour  of  results  using  a  discre¬ 
tization  of  150  X  3  rings,  and  varying  the  length  of  the  tube. 


Fig.  5  Induced  voltage  vs.  coil  velocity,  tube  witliout  defect 


Fig.  5  shows  the  calculated  and  analytical  voltage  on  the 
pick-up  coil  for  different  coil  velocities  and  three  different 
tube  lengths.  The  results  show  that  the  length  of  the  tube  has 


Fig.  7  Induced  current  modulus  in  the  tube  vs.  distance  from  tire 
exciting  coil  on  three  radii,  for  v  =  0. 


Fig.  8  Induced  current  modulus  in  the  tube  vs.  distance  from  tire 
exciting  coil  on  three  radii,  for  v  =  100  m/sec. 

The  effects  of  velocity  are  clearly  shown  by  comparison 
of  the  analytically  and  numerically  calculated  induced  cur¬ 
rents  in  the  tube.  Fig.  7  and  8  show  the  induced  currents  at 
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three  different  radii  in  the  tube  for  v=0  and  v=100  ni/sec  for 
the  previous  geometry  and  discretization  level. 

Next,  we  study  the  geometry  shown  in  fig.  9  considering 
a  conductivity  a=  2-10^  (S/m)  and  a  constant  relative  per¬ 
meability  Pr=1000,  the  e.\citing  signal  frequency  is  50  Hz. 
For  configurations  which  includes  an  axisymmetric  defect, 
there  is  an  analytical  solution  only  for  a  non  ferromagnetic 
tube  [5],  or  for  a  ferromagnetic  tube  at  zero  coils  velocity  [7], 


Fig.  9  RFEC  geometry  with  an  axisymmetric  defects. 

Fig.  10  and  1 1  show  good  agreement  between  the  com¬ 
puted  and  analytical  induced  voltage  and  phase  for  a  defect 
size  of  2  X  4  mm  and  different  relative  positions  between  the 
defect  and  the  coils.  In  the  following  figures  distance  zero 
mpqns  that  the  pick  up  coil  is  exactly  under  the  defect. 


Fig.  10  Induced  voltage  vs  coil  position,  tube  witli  defect 


Fig.  1 1  Phase  of  the  induced  voltage  vs.  coil  position,  tube  with 
defect 


As  a  conclusion  from  the  tests  above,  we  observe  that  the 
characteristic  "double  bump"  signal  was  only  clearly  detect¬ 
able  if  the  tube  material  has  a  high  conductivity;  copper  for 
instance.  The  influence  of  tube  length  and  discretization 
level  on  the  results  calculated  for  the  cases  shovm  in  figures 
7,  8.  10  and  11  is  not  shown  here  but  these  show  the  same 
trends  as  the  results  in  fig.  4  and  5.  The  influence  of  the 
tube  length  is  more  pronounced  than  that  due  to  discretiza¬ 
tion.  This  means  that  the  proposed  method  does  not  require 
fine  discretization  for  good  results,  even  in  presence  of 
defects.  One  purpose  of  this  method  was  to  show  that  the 
RFEC  is  a  direct  induction  effect.  To  do  so  we  calculated  the 
induced  voltage  for  tubes  with  wall  thickness  of  7  cm.  and 
12  cm. 


Fig.  12  Induced  voltage  vs  coil  position,  thick  tube  with  inner 
defect 


Figure  12  shows  that  the  induced  voltage  is  still  sensitive 
to  the  presence  of  an  internal  defect  as  the  pick-up  coil 
passes  under  the  defect.  In  addition,  the  magnitude  of  the 
induced  voltage  is  nearly  equal  for  the  two  thicknesses 
tested,  indicating  that  sensitivity  is  unrelated  to  wall  thick¬ 
ness.  Figures  13,  and  14  show  the  transient  induced  voltage 
for  different  frequencies  and  probe  velocity. 


Fig.  13  hiduced  voltage  vs.  coil  position,  f=50  Hz,  v=5ra/sec 


distance  [ml  _  distance  [m] 

Fig.  14  Muced  voltage  vs.  coil  position,  f=50  Hz,  v=10ni/sec  Fig.  16  Induced  voltage  vs.  coil  position,  f=7850  Hz,  v=50m/sec 


The  induced  voltage  has  a  typical  transient  behaviour, 
even  for  the  lower  velocities  considered  here.  This  behaviour 
cannot  be  related  to  steady  state  parameters.  Furthermore, 
there  is  also  a  dependence  between  coil  velocity  and  e.\cita- 
tion  frequency  and  detectability  of  the  defect.  In  fact,  at  a 
frequency  of  50  Hz,  and  a  velocity  of  10  m/sec,  the  pick-up 
signal  is  distorted,  as  shown  in  fig.  12,  while  at  a  velocity  5 
m/sec  the  amplitude  variation  is  clearer,  and  the  defect  can 
be  discerned.  Therefore,  the  conclusion  is  that  there  must  be 
a  high  ffequency/velocity  ratio  for  clear  defect  detection. 
This  is  also  pointed  out  in  the  figures  15  and  16. 


distance  (m} 

Fig.  15  Induced  voltage  vs.  coil  position,  f=785  Hz,  v=5m/sec 


Figure  15  shows  a  clearer  variation  of  the  pick-up  signal 
in  comparison  with  Fig.  12.  This  is  accomplished  by  increas¬ 
ing  the  frequency  from  50  Hz  to  785  Hz.  Fig.  16  shows  that 
the  sensitivity  is  unchanged  if  the  frequency/velocily  ratio  is 
kept  constant  but  with  a  frequency  and  velocity  10  times 
higher. 


4.  CONCLUSIONS 

The  integral  formulation  presented  here  allows  simple 
modeling  of  RFEC  tests,  including  defects.  Computation 
time  is  relatively  short  even  in  the  presence  of  moving  ob¬ 
jects.  The  method  has  been  validated  first  by  comparison 
with  analytical  results.  Application  of  the  formulation  for 
the  analysis  of  the  response  of  a  RFEC  testing  apparatus  has 
highlighted  the  importance  of  several  aspects  of  the  tran¬ 
sient  phenomena. 

In  particular,  the  output  signal  must  be  transient  for 
complete  understanding  of  defect  detection  process.  Also, 
the  importance  of  the  relation  between  frequency  of  the  ex¬ 
citing  signal  and  coil  velocity  was  pointed  out. 
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(  )  $116 

EUROPE,  FORMER  USSR 
TURKEY,  SCANDINAVIA 

(  )  $68 

(  )  $78 

(  )  $116 

ASIA,  AFRICA,  MIDDLE 

EAST  &  PACIFIC  RIM 

(  )  $68 

(  )  $85 

(  )  $116 

FULL-TIME  STUDENT  RATE  IS  $25  FOR  ALL  COUNTRIES 


♦THIS  SECTION  IS  TO  BE  USED  ONLY  IF  PAYING  BY  CREDIT  CARD  &  CARD  IS  NOT  YOUR  OWN. 
IT  IS  NOW  MANDATORY  THAT  THE  CARD  HOLDER  PRINTS  AND  SIGN  HIS/HER  NAME  AND 

LIST  HIS/HER  COMPLETE  ADDRESS* 


PRINT  FIRST  AND  LAST  NABIE  OF  CARD  HOLDER  SIGNATURE  OF  CARD  HOLDER 


BiAIUNG  ADDRESS 

BIAILING  ADDRESS(cont) 


METHOD  OF  PAYMENT 

□ 

A  bank  check  for  the  total  eimount  is  enclosed. 

□ 

Traveler’s  checks  for  the  total  amount  are  enclosed.® 

□ 

International  Money  Order  is  enclosed.® 

□ 

Charge  to:  □  MasterCard  nVisa.  n  Discover  O  Amex.® 

Card 

No. 


Card  Exp.  Date: 
Mo. _  Year 


MAKE  PAYABLE  TO  "ACES"  and  send  to:  RICHARD  W.  ADLER,  EXEC.  OFFICER, 

NAVAL  POSTGRADUATE  SCHOOL,  ECE  DEPT.,  CODE  EC/AB,  833  DYER  ROAD,  ROOM  437,  MONTEREY,  CA  93943-5121 

Non-USA  participants  may  remit  via  (1)  Bank  Checks,  if  (a)  drawn  on  a  U.S.  Bank,  (b)  have  bank  address,  (c)  contain 
series  of  (9)  digit  mandatory  routing  numbers;  (2) Traveler’s  Checks  (in  U.S.  $$) ;  (3)  International  Money  Order  drawn 
in  U.S.  funds,  payable  in  U.S.;  (4)  Credit  Cards:  Visa,  Master  Card,  Discover  Card,  Amex. 

Total  Remittance  (U.S.  Dollars  Only)  $ _  July  1996 
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ADVERTISING  RATES 


FEE 

PRINTED  SIZE 

Full  page 

$200. 

7.5"  X  10.0” 

1/2  page 

$100. 

7.5"  X  4.7"  or 

3.5"  X  10.0” 

1/4  page 

$  50 

3  5"  X  4.7" 

All  ads  must  be  camera  ready  copy. 

Ad  deadlines  are  same  as  Newsletter  copy  deadlines. 

Place  ads  with  Ray  Perez,  Newsletter  Editor,  Martin  Marietta  Astronautics, 

MS  58700,  PO  Box  179,  Denver,  CO  80201,  USA. 
reject  ads. 

The  editor  reserves  the  right  to 

ACES  NEWSLETTER  AND  JOURNAL  COPY  INFORMATION 


Issue  Copy  Deadline 

March  Januaiy  13 

July  May  25 

November  September  25 
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INFORMATION  FOR  AUTHORS 
PUBLICATION  CRITERIA 


Each  paper  is  required  to  manifest  some  relation  to  applied 
computational  electromagnetics.  Papers  may  address 
general  issues  in  applied  computational  electromagnet¬ 
ics,  or  they  may  focus  on  specific  applications,  tech¬ 
niques,  codes,  or  computational  issues.  While  the 
following  list  is  not  exhaustive,  each  paper  will  generally 
relate  to  at  least  one  of  these  areas: 

1.  Code  validation.  This  is  done  using  internal  checks  or 
experimental,  analytical  or  other  computational  data. 
Measured  data  of  potential  utility  to  code  validation  efforts 
will  also  be  considered  for  publication. 

2.  Code  performance  analysis.  This  usually  involves 
identification  of  numerical  accuracy  or  other  limitations, 
solution  convergence,  numerical  and  physical  modeling 
error,  and  parameter  tradeoffs.  However,  it  is  also 
permissible  to  address  issues  such  as  ease-of-use,  set-up 
time,  run  time,  special  outputs,  or  other  special  features. 

3.  Computational  studies  of  basic  physics.  This  involves 
using  a  code,  algorithm,  or  computational  technique  to 
simulate  reality  in  such  a  way  that  better  or  new  physical 
insight  or  understanding  is  achieved. 

4.  New  computational  techniques,  or  new  applications  for 
existing  computational  techniques  or  codes. 

5.  "Tricks  of  the  trade"  in  selecting  and  applying  codes 
and  techniques. 

6.  New  codes,  algorithms,  code  enhancement,  and  code 
fixes.  This  category  is  self-explanatory  but  includes 
significant  changes  to  existing  codes,  such  as  applicability 
extensions,  algorithm  optimization,  problem  correction, 
limitation  removal,  or  other  performance  improvement. 
Note:  Code  (or  algorithm)  capability  descriptions  are 
not  acceptable,  unless  they  contain  sufficient  technical 
material  to  justify  consideration. 

7.  Code  input/output  issues.  This  normally  involves 
innovations  in  input  (such  as  input  geometry 
standardization,  automatic  mesh  generation,  or  computer- 
aided  design)  or  in  output  (whether  it  be  tabular,  graphical, 
statistical,  Fourier-transformed,  or  otherwise  signal- 
processed).  Material  dealing  with  input/output  database 
management,  output  interpretation,  or  other  input/output 
issues  will  also  be  considered  for  publication. 

8.  Computer  hardware  issues.  This  is  the  category  for 
analysis  of  hardware  capabilities  and  limitations  in  meeting 


various  types  of  electromagnetics  computational  require¬ 
ments.  Vector  and  parallel  computational  techniques  and 
implementation  are  of  particular  interest. 

Applications  of  interest  include,  but  are  not  limited  to, 
antennas  (and  their  electromagnetic  environments), 
networks,  static  fields,  radar  cross  section,  shielding, 
radiation  hazards,  biological  effects,  electromagnetic  pulse 
(EMP),  electromagnetic  interference  (EMI),  electromagnet¬ 
ic  compatibility  (EMC),  power  transmission,  charge 
transport,  dielectric  and  magnetic  materials,  microwave 
components,  MMIC  technology,  remote  sensing  and  geo¬ 
physics,  communications  systems,  fiber  optics,  plasmas, 
particle  accelerators,  generators  and  motors,  electromagnet¬ 
ic  wave  propagation,  non-destructive  evaluation,  eddy 
currents,  and  inverse  scattering. 

Techniques  of  interest  include  frequency-domain  and 
time-domain  techniques,  integral  equation  and  differential 
equation  techniques,  diffraction  theories,  physical  optics, 
moment  methods,  finite  differences  and  finite  element 
techniques,  modal  expansions,  perturbation  methods,  and 
hybrid  methods.  This  list  is  not  exhaustive. 

A  unique  feature  of  the  Journal  is  the  publication  of 
unsuccessful  efforts  in  applied  computational 
electromagnetics.  Publication  of  such  material  provides  a 
means  to  discuss  problem  areas  in  electromagnetic  model¬ 
ing.  Material  representing  an  unsuccessful  application  or 
negative  results  in  computational  electromagnetics  will  be 
considered  for  publication  only  if  a  reasonable  expectation 
of  success  (and  a  reasonable  effort)  are  reflected. 
Moreover,  such  material  must  represent  a  problem  area  of 
potential  interest  to  the  ACES  membership. 

Where  possible  and  appropriate,  authors  are  required  to 
provide  statements  of  quantitative  accuracy  for  measured 
and/or  computed  data.  This  issue  is  discussed  in  "Accuracy 
&  Publication:  Requiring  quantitative  accuracy  statements 
to  accompany  data",  by  E.K.  Miller,  ACES  Newsletter,  Vol. 
9,  No.  3,  pp.  23-29,  1994,  ISBN  1056-9170. 

EDITORIAL  REVIEW 

In  order  to  ensure  an  appropriate  level  of  quality  control, 
papers  are  refereed.  They  are  reviewed  both  for  technical 
correctness  and  for  adherence  to  the  listed  guidelines 
regarding  information  content.  Authors  should  submit  the 
initial  manuscript  in  draft  form  so  that  any  suggested 
changes  can  be  made  before  the  photo-ready  copy  is 
prepared  for  publication. 
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INFORMATION  FOR  AUTHORS 
STYLE  FOR  CAMERA  READY  COPY 


The  ACES  Journal  is  flexible,  within  reason,  in  regard  to 
style.  However,  certain  requirements  are  in  effect: 

1.  The  paper  title  should  NOT  be  placed  on  a  separate 
page.  The  title,  author(s),  abstract,  and  (space  permitting) 
beginning  of  the  paper  itself  should  all  be  on  the  first  page. 
The  title,  author(s),  and  author  affiliations  should  be 
centered  (center-justified)  on  the  first  page. 

2.  An  abstract  is  REQUIRED.  The  abstract  should  state 
the  computer  codes,  computational  techniques,  and 
applications  discussed  in  the  paper  (as  applicable)  and 
should  otherwise  be  usable  by  technical  abstracting  and 
indexing  services. 

3.  Either  British  English  or  American  English  spellings 
may  be  used,  provided  that  each  word  is  spelled  consistently 
throughout  the  paper. 

4.  Any  commonly-accepted  format  for  referencing  is 
permitted,  provided  that  internal  consistency  of  format  is 
maintained.  As  a  guideline  for  authors  who  have  no  other 
preference,  we  recommend  that  references  be  given  by 
author(s)  name  and  year  in  the  body  of  the  paper  (with 
alphabetical  listing  of  all  references  at  the  end  of  the  paper). 
Titles  of  Journals,  monographs,  and  similar  publications 
should  be  in  boldface  or  italic  font  or  should  be  underlined. 
Titles  of  papers  or  articles  should  be  in  quotation  marks. 

5.  Internal  consistency  shall  also  be  maintained  for  other 
elements  of  style,  such  as  equation  numbering.  As  a 
guideline  for  authors  who  have  no  other  preference,  we 
suggest  that  equation  numbers  be  placed  in  parentheses  at 
the  right  column  margin. 

6.  The  intent  and  ^meaning  of  all  text  must  be  clear.  For 
authors  who  are  NOT  masters  of  the  English  language,  the 
ACES  Editorial  Staff  will  provide  assistance  with  grammar 
(subject  to  clarity  of  intent  and  meaning). 

7.  Unused  space  should  be  minimized.  Sections  and 
subsections  should  not  normally  begin  on  a  new  page. 

MATERIAL,  SUBMITTAL  FORMAT  AND 
PROCEDURE 

The  preferred  format  for  submission  and  subsequent  review, 
is  12  point  font  or  12  cpi,  double  line  spacing  and  single 
column  per  page.  Four  copies  of  all  submissions  should  be 
sent  to  the  Editor-in-Chief.  Each  submission  must  be 
accompanied  by  a  covering  letter.  The  letter  should  include 
the  name,  address,  and  telephone  and/or  fax  number  and/or 
e-mail  address  of  at  least  one  of  the  authors. 


Only  camera-ready  original  copies  are  accepted  for 
publication.  The  term  "camera'- ready”  means  that  the 
material  is  neat,  legible,  and  reprodedble.  The  preferred 
font  style  is  Times  Roman  10  point  (or  equivalent)  such  as 
that  used  in  this  text.  A  double  column  format  similar  to 
that  used  here  is  preferred.  No  author's  work  will  be 
tnrned  down  once  it  has  been  accepted  because  of  an 
inability  to  meet  the  requirements  concerning  fonts  and 
formats  Full  details  are  sent  to  the  author(s)  with  the  letter 
of  acceptance. 

There  is  NO  requirement  for  India  ink  or  for  special  paper; 
any  plain  white  paper  may  be  used.  However,  faded  lines  on 
figures  and  white  streaks  along  fold  lines  should  be  avoided. 
Original  figures  -  even  paste-ups  -  are  preferred  over 
"nth-generation"  photocopies.  These  original  figures  will 
be  returned  if  you  so  request. 

While  ACES  reserves  the  right  to  re-type  any  submitted 
material,  this  is  not  generally  done. 

PUBLICATION  CHARGES 

ACES  members  are  allowed  12  pages  per  paper  without 
charge;  non-members  are  allowed  8  pages  per  paper  without 
charge.  Mandatory  page  charges  of  $75  a  page  apply  to  all 
pages  in  excess  of  12  for  members  or  8  for  non-members. 
Voluntary  page  charges  are  requested  for  the  free  (12  or  8) 
pages,  but  are  NOT  mandatory  or  required  for  publication. 
A  priority  courtesy  guideline,  which  favors  members, 
applies  to  paper  backlogs.  Full  details  are  available  from  the 
Editor-in-Chief. 

COPYRIGHTS  AND  RELEASES 

Each  primary  author  must  sign  a  copyright  form  and  obtain 
a  release  from  his/her  organization  vesting  the  copyright 
with  ACES.  Forms  will  be  provided  by  ACES.  Both  the 
author  and  his/her  organization  are  allowed  to  use  the 
copyrighted  material  freely  for  their  own  private  purposes. 

Permission  is  granted  to  quote  short  passages  and  reproduce 
figures  and  tables  from  an  ACES  Journal  issue  provided  the 
source  is  cited.  Copies  of  ACES  Journal  articles  may  be 
made  in  accordance  with  usage  permitted  by  Sections  107 
or  108  of  the  U.S.  Copyright  Law.  This  consent  does  not 
extend  to  other  kinds  of  copying,  such  as  for  general 
distribution,  for  advertising  or  promotional  purposes,  for 
creating  new  collective  works,  or  for  resale.  The 
reproduction  of  multiple  copies  and  the  use  of  articles  or 
extracts  for  commercial  purposes  require  the  consent  of  the 
author  and  specific  permission  from  ACES.  Institutional 
members  are  allowed  to  copy  any  ACES  Journal  issue  for 
their  internal  distribution  only. 


